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)nt1sp1ece: A on form flow (left) beneath a cylindrical spreader loses 
stability (center: onset of Instability ) to a transversely 
period c ribbing, disturbance (right) as we sjee ?S!S!« 
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CHAPTER 1 - PREL UDE 

i 

1.1. Introduction 

Coating the surface of a solid or flexible sheet uniformly with a 
thin film of liquid Is Important In a variety of scientific and indus- 
trial applications- Imperfections In the coating, especially as the 
soeed of coating, is increased or the coating thickness is decreased, are 
mostly due to now instabilities and various wetting and edge effects. 
The goal of this thesis fs to understand, by means of theoretical analy- 
sis and computer simulation,, some .of the Instabilities in thin liquid 
films as they are deposited, and the nonunl f orraltles of film thicKness 
that result. These are matters of hydrodynamic instability and bifur- 
cation or branching of families of flow states, for which conventional 
theoretical approaches ».nd more novel computer-based methods are 
advancing. At the center of this thesis -f-s the extension to film-flow 
problems of Galerk1n/f1n1te element nsthods of stability analysis 
recently developed by Brown and Scrlven (1980a) and also reported in 
Brown's thesis (1979) for liquid drops 1n hydrostatic equilibrium, which 
are simpler to treat than flowing liquid films. 

This thesis is built upon a foundation laid by a host of contribu- 
tors. As mentioned 1n the foregoing paragraph, one of the most Important 
of them was Brown. Much of the success of this thesis is owed to the 
experience and insights gained by Orr (1976), Sllliman (1979, 1980 with 
Scrlven), and Saito (1981 with Scrlven) by their efforts to develop and 
improve Galerkin/finlte element algorithms for solving steady, two- 
dimensional, free surface How problems. These earlier works ire sum- 



2 

marl2fcu In a recent article by Mstler and Scrlven U982), who also 
describe several Improvements and generalizations that have been recently 
d*i?1sed. Finally, an asymptotic analysis worked out by H1gg1ns (1980, 
1=32) showing how film flows approach a final regime Is the cornerstone 
of Chapter 2 and plays an Important role In the methcfds described 1n 
Chapters 4 and 5. 

Understanding the mechanisms of instability of uniform flowing 
films 1s not only a scientific challenge o.; the frontier zf capillary 
hydrodynamics: such understanding will point the way to avoiding or at ■ 
least controlling imperfections a'nd incipient Film breakdowns wherever 
they are not absolutely inherent in a ccating system. 

Capillary H ydrodynamics and Coating Rows . By their very nature, 
coating flows are free surface Hows. Free surfaces are menisci, i.e. 
liquid/gas interfaces, and their presence greatly complicates flow field 
analyses because their location, and thus the flow domain, 1s not known 
a priori. Menisci are In reality layers of finite thickness in which 
density varies rapidly, from that of bulk liquid to that of bulk gas. 
However, the thicknesses are usually no more than tens of Angstroms and 
therefore are nearly Infinitesimal In comparison with fluid mechanical 
length scales. It is convenient to approximate menisci as satheaatical 
surfaces. The steep density gradients at a meniscus give rise to a 
property known as surface tension which produces an isotropic capillary 
stress within the meniscus not unlike that In a stretched rubber nem- 
brane. 

Another feature comnon to. coating flows 1s that they contain three- 



phase contact lines where menisci Intersect solid surfaces. Contact 
* lines come <n two varieties: (U they are static If their position 
relative to the solid boundary is steady; (2) they are dynamic if they 
translate relative to the solid boundary. Oynamic contact lines are 
also called wetting lines. Every coating flow has a wetting line 
because liquid must come In contact with a dry web, I.e. the sheet being 
coated. Most coating Hows also have static contact lines at the place 
where the liquid being coated leaves the coating die. 

Figure 1.1 portrays five commonly used devices; apart from the 
knife coater they are suitable for premetered, precision coating. The 
slot coater, also called the extrusion coater, consists of a solid die 
through which liquid is extruded at a metered rate. The die may be 
parallel or inclined to the moving web that is supported on the opposite 
side by a back-up roll, or in some cases stretched between two rolls not 
far apart. The downstream separating contact line can be pinned at a 
geometrical or compositional discontinuity, *$ it Is 1n Figure 1.1, or 
1t can be free to move along the solid die. The knife coater is essen- 
tially the *>wnstre«ia portion of the slot coater. The roll coater con- 
sists of two rolls, one of which usually supports a web. The liquid 
entrained between the rolls splits at the meniscus so that part adheres 
to the web and the remainder Is carried away by one of the rolls. A 
cannon variation on roll coating Is reverse roll coating in which the 
two rolls are made to corotate instead of counterrotate. Coyle et al. 
{1332a, b) have Investigated roll coating In both modes of operation by 
finite element analysis. The slide coater is made up of a die which Is 



inclined with respect to gravity and a web backed by u rcll which Is 
separated from the die by a small gap. Liquid Hows down th* die, or 
slide, by gravity and Is transferred to the web. 

Each of the coaters pictured In Figure l.l has finite width, but 
because the rtdth'ls usually ouch greater than the other length scales 
of the coater It Is coureon to make the mathematical Idealization that 
coating Hows are transversely uniform. I.e. two T d1mensional . indeed, 
one of the Important tasks of coating engineers Is to design coaters so 
that edge effects are minimized in order to make the wliole width of the 
coating as uniform as possible. Ideally, coating flows are not only 
transversely uniform but also steady, for otherwise the coating thickness 
would not be uniform In the strearawlse direction. However. 1n practice 
it Is never possible to maintain complete transverse or streamwlse uni- 
formity. The degree to which uniformity can be achieved Is controlled 
not only by the amplitude and frequency of disturbances such as mechani- 
cal vibrations and pressure fluctuations, but also by the sensitivity of 
the flow to these disturbances or even to infinitesimal disturbances due 
to thermal fluctuations. 

1.2. Governing Equations and Boundary Conditions 

Coating Hows obey well-known equations. which are devised from phys- 
ical principles, inese principles are conservation of momentum and mass, 
and they are expressed mathematically for Incompressible Hows by the 
Navler-Stokes equations: 
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Re[3u/3t + 7.(uu)] * 7-T + St (1.2.1) 

0aV 'i • (1.2.2) 

Her- s!:: vector representing local velocity, T 1s a dyadic or tensor 
expressing the local state of stress, and e^ 1s a unit vector pointing 
In the direction of gravity. Reynolds number, Re = Ud/v. measures the 
relative Importance of Inertia to viscous force; Stokes number, St = 
gd M measures the relative importance of the force due to gravitational 
acceleration to viscous force. Here U Is a characteristic velocity, 
e.g. web speed; d Is a characteristic length, e.g. gap width; v Is kine- 
matic viscosity; and g Is local gravitational acceleration. This thesis 
deals exclusively with flows in which the liquids are Newtonian and have 
unlfor* properties. In this case T is given by: 

T » Tu + (7u) T - pi 

The superscript T Indicates the transpose of the dyadic or tensor, and 
P Is the Isotropic portion of the state of stress. The tensor I = ij_ ♦ 
JJ + Jclc, and 1_, J, and £ are unit vectors denoting the x-, y-, and z- 
dlrectlons in a Cartesian coordinate system. 

Ordinarily liquids adhere to. solid walls, or very nearly so. When 
this Is so, the liquid velocity at a wall Is identical to wall velocity 
U:' 

U a U ' , 

- - (1.2.3) 
This adherence condition appears to be violated, however, in some cases 



where shear stress becomes very high at the wad, such as at contact 
lines and especially at dynamic contact lines as hypothesized by Huh. and 
Scriven (1971) and observed experimentally by Kraynik and Schowalter 
(1977) and others.' In a forerunner of this thesfs, Slll.iraan (1979, 1930 
with Scriven) showed that allowing some slip at vettfmj'Mnes relieves 
the nonlntegrable mathematical singularities which would otherwise be 
present. One way to allow jlip is to replace (1.2.3) by 

u * f(s) u^ + (.1 - f( s )] U . (1.2.4) 

Here ^ Is the velocity of the contact line and f(s) is a function of 
arclength s measured from the contact line along the wall so that f(0) 
3 1. He) = O f and f(s) * 0 when s > c for some small distance e. Equa- 
tion (1.2.4) serves the same purpose as the slip boundary condition pro- 
posed by Navler (1827), but has the advantage of being easier to employ. 
However, (1.2.41 has the disadvantage that it can lead to aphyslcal 
results - I.e. 1t can predict momentum transport up a velocity discon- 
tinuity, 1n violation of the second law of thermodynamics - unless one 
is careful to select a function f(s) which decreases. monotonically on 
the Interval [0,1]. 

Three boundary conditions apply at a free surface. The first two 
are derived from a force balance there; the third is a statement of mass 
conservation. The difference In the total normal stress exerted on an 
interface by the liquid phase (L) and gas phase (G) roust be equal to the 
capillary stress: 



0 * (T L - T G ):nn - Ca~ l 2H 



(1.2.5) 
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Here n Is the outward unit normal to the interface, H Is the mean cur- 
vature of the interface, and Ca s ull/c is the capillary number, <x and u 
being respectively the surface tension and 4/namic viscosity. A useful 
approximation is that the gas is invlsdd, which is usually valid when* 
ever the liquid viscosity is much higher than that of the gas* If in 
addition the gas pressure is taken as a datum, f vanishes altogether. 
When surface tension Is uniform, as it is taken to be throughout this 
thesis, the corwUion that shear stresses balance at an interface leads 
to the condition 

0 = [fT L - fj.nl xn ■ (1 . 2i6 , 

Finally, the condition that the -interface be a material surface is the 
kinematic equation 

* 7( *- h) (1.2-7) 
Here y * h represents the location of the free surface. 

1.3. Steady Two-Olmensional Ccatlng Flows 

The starting point in aty analysis of a coating flow Is to look for 
steady solutions of the governing equations and boundary conditions as 
given in {1.2.lj-(1.2.7). In most. instances transverse uniformity of the 
coating Is sought and so Hows are idealized as being two-dimensional. 
Even with these simplifications, constructing solutions by conventional 
mathematical methods Is not a trivial matter unless further approxima- 
tions are made. 
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A set of approximations, collectively termed the lubrication approx- 
imation, can be made In order te simplify the governing equations* so that 
they can be handled by conventional mathematical methods. The approxi- 
mations are: (15 Reynolds number times a characteristic slope of the die 
is small, I.e.. web speed times gap thickness times a characteristic slope 
divided by kinematic viscosity Is small; (2) flow is nearly rectilinear, 
i.e. velocity profiles are superpositions of shear-driven Couette and 
pressure-driven Pbiseullle components, because the flow is confined to a 
channel which 1s nearly parallel; (3) liquid properties are uniform; and 
(4) acceleration due to gravity 1s' unimportant (removable restriction). 
The lubrication approximation Is used in Chapter 3 to investigate slot 
or knife coating flows when the die has an exponential shape. 

The chief drawback of lubrication theory is that the approximation 
breaks down In the region near a separating contact line because the 
assumption of nearly rectilinear flow falfs there. The usual remedy is 
to represent toe flow field only as far downstream as the meniscus, and 
to impose a set of ad hoc conditions there. However, it Is not certain 
which Is worse: the remedy or the disease. The problem is that a 
collection of ad hoc boundary conditions have been postulated, and they 
each predict a different meniscus location as shown in Figure 1.2. 

The alternative is to construct solutions to Equations (1.2.1)- 
(1.2.7) by the-Galerkin/flnite element method. The i ox* uf finite ele- 
ment analysis of steady flow states Is to repres- the velocity and 
pressure fields and the free surface location i --rws of particularly 
convenient sets of basis functions. These are each nonzero on only a 
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subdorsal n of the flow domain, In contrast to the traditional functions 
of mathematical physics, which span ait entire domain. Moreover, the basis 
function on ea<* s'?bd«mln is a simple polynomial, so that more and more 
accurate estimates of the flow field and free surface shape can be cal- 
culated by dividing the domain Into more andnore subdomains. However, 
as the number increases, more computations are required to find the 
coefficients In all the polynomial basis functions. so that (1) mass and 
momentum are conserved and boundary conditions are satisfied on the 
scale of the subdomains, and (2) the estimate of the solution has the 
desired degree of continuity. Mot only is the Idea straightforward; it 
is also consistent with rigorous mathematics. But its implementation Is 
so tedious as to be Impracticable without sophisticated matrix calcula- 
tions for the large systems of algebraic equations to which the finite 
element method ultimately leads. Such calculations have become feasible 
only with the advent of large, fast digital computers. Analyses of 
coating flow problems by the finite element method are expensive. Thus, 
It Is imperative to find efficient techniques arid to look for ways of 
reducing the number of equations and unknowns without sacrificing 
accuracy. 

One way to reduce the number of equations and unknowns without sac- 
rificing accuracy Is to impose conditions at inflow and outflow bound- 
aries which represent well the physics there so that the computational 
domain can be shrunk. The idea Is to split up the domain Into three 
rones, as illustrated in Figure 1.3. In the upstream channel zone the 
flow Is near]y parabolic so that perturbation analysis can be used to 
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solve the governing equations there. Likewise, flow is nearly solid 
body translation In the downstream free surface zone so that perturbation 
analysis car. be used there as well. These analyses were carried out by 
Wilson (1969) and Htggins (1980, 1982) for the channel and. free surface 
zones, respectively. The Galerkin/finlte element method Is used to rep- 
resent the flow in the forming zone. Matching conditions between forming 
zone and channel zone, I.e. at the Inflow boundary of the forming zone, 
and between forming zone and free surface zone, i.e. at the outflow 
boundary of the forming zone, take the form of Robin conditions. I.e. 
boundary conditions of the third kind. The details are set out 1n 
Chapter 4. The main findings there are that at the inflow boundary the 
Robin condition offers' little advantage In ccsparl son. with more conven- 
tional Olrichlet or Heumann conditions, i.e. boundary conditions In 
Mch either the dependent variables or derivatives of the dependent 
variables are specified. On the other hand~ the length of the forming 
zone measured downstream of the contact line can be reduced by as much 
as fourfold without sacrificing computational accuracy by using a Robin 
condition there instead of a Oirlchlet or Neumann condition, as shown in 
Figure 1.4. 

1.4. Steady Three-Dlmensional Film Hows 

Because coatings are never flawless In practice, coating engineers 
are forced to alii only for uniformity within some tolerance. Surface 
leveling can play an Important role by reducing the amplitude of coating 
imperfections. Estimates of leveling rates could aid significantly in 
the design of coating dies by making clear which disturbances level so 
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slowly that they cannot be tolerated. 

In Chapter 2 biggins 1 (1980. 1SS2) asymptotic analysis of the two-, 
dimensional development of vlscocaplllary film flows 1s extended to 
include three-dimensional disturbances. Three-dimensional disturbances 
are representable In a doubly infinite' expansion in domain-spanning basis 
functions, which are the eigenvectors of the elgenproblem to which the 
perturbation analysis reduces. However, one mode generally persists 
further downstream than any of the others, especially when the upstream 
disturbance consists of a single Fourier mode, so that only this dominant 
mode need be considered If one's primary Interest is in the flow far 
downstream. 

The main results are the dominant decay rates as functions of 
Reynolds number, capillary number, and wavenumber, I.e. the inverse of 
the transverse wavelength of the disturbance, as contained In Tables 
2.1-2.5. An Important distinction with the earlier work of Hlggins 
(1980, 1982) is that the dominant mode can be one of damped standing 
waves when disturbances are three-.dlmenslonal , as Illustrated in Figure 
1.5a, whereas Hlggins found that the dominant mode always decays expo- 
nentially when disturbances are two-dimensional. Figure 1.5b shows an 
example of a three-dimensional disturbance that decays exponentially. 

An important use of the results of Chapter 2 Is to extend the Robin 
condition tested in Chapter 4 to a three-dimensional analogue: 



Here T Is the three-dimensional state of stress tensor, u is the three- 
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Figure 1.5 Examples of upstream disturbances that decay with downstream 
length (a) to form damped standing waves, and (b) exponen- 
tially. Damped standing waves only appear when the distur- 
bance 1s three-dimensional. Both cases are at Reynolds 
number Re = 100. and wavenumber N = 1. Case (a) is at 
capillary number Ca = 100. and case (b) is at capillary 
number Ca = 1 . 
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dimensional velocity vector, and A Is a proportionality tensor that is 
constructed from the asymptotic solution. Equation (1.4.1) Is espe- 
cially useful in stability analyses of steady two-dimensional Hows to 
three-dimensional disturbances. 

1.5. Stability of steady Two-Olmensional Coating Flows 

Even when coating die Imperfections *r misalignment, mechanical 
vibrations, and pressure fluctuations are held In check, there is no 
guarantee that a coating win be uniform. For if the flow Is fluld- 
mechanically unstable, even infinitesimal disturbances due to thermal 
fluctuations or the like can be amplified to finite size. It is there- 
fore Insufficient to show that steady two-dimensional solutions of the 
governing equations exist; It Is also necessary to show that they are 
stable . 

One of the most prevalent of coating now Instabilities 1s the rib- 
binj Instability. Ribbing is a steady pattern of waves which leads to 
nonunlformity of the coating in the transverse direction. Each of the 
coaters shown In Figure 1.1 1s known to produce ribbed coatings under 
certain conditions. Usually, the faster the coating speed, the lower 
the surface tension, the higher the liquid viscosity, or the thinner the 
coating, the less stable is a transversely uniform flow to a ribbing 
disturbance. 

Figure 1.6 Is a photograph taken by Dennis Coyle showing ribbing 
In a two-roll setup similar to the one shown In Figure 1.1c except that 
the two rolls were ha If- lomersed In a bath. Even though the apparatus 



Figure 1.6. The ribbing instability in a two-roll set-i 
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is nearly transversely symmetric (aside from effects of the edges), the 
liquid film that adheres to the rolls Is far from being uniform. Instead, 
the film thickness Is remarkably periodic in the transverse, or axial, 
direction. 

Figure 1.7 shows, sone results of a visualization experiment that was 
performed at a research facility of Minnesota Mining and Manufacturing 
with the aid of Gary Zvan. The coater that was used Is equivalent to 
the knife- coater. shown 1n Figure 1.1 except that the upper die was a cyl- 
inder of 3.5 csntiaeter diaueter. Th a minimum gap between the cylinder 
and the web was set to 0.025 centimeters. The liquid used was glycerin, 
which had a measured viscosity of 1250 centl poise and surface tension of 
69.5 dynes per centimeter. Coatings were made by depositing a small 
puddle of liquid upstream of the cylinder so that liquid was spread by 
the cylinder onto the traveling web. The coating was then photographed 
by a camera mounted directly above the web. A sheet of graph paper 
having millimeter rulings was held slightly below and parallel to the 
transparent web so that nonunlformitles could be more readily observed. 
Web motion Is from top to bottom In the three frames In Figure 1.7. The 
coating appeared uniform at the slowest web speed shown In Figure 1.7a. 
At the intermediate speed, a pattern of ribs similar to those In Figure 
1.6 was observed, although they are barely noticeable In Figure 1.7b. 
Figure 1.7c shows a distinct pattern of ribs which have a dimenslonless 
wavenumber of about 0.53 (made dimenslonless by the minimum gap clear- 
ance), i.e. the dimensional wavenumber is 2e/wavelength. 
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Stability of Coatin g Rows Under the Lubrication Approximation . 
The stability of a steady two-dimensional coating flow can be estimated 
by using the lubrication approximation to account not only for the steady 
two-dimensional flow but also for possible unsteady three-dimensional 
disturbances to It.' Pearson (I960) was the first to show how this can 
be done in his classical analysis of ribbing in a wedge-shaped spreader. 
However, the same question arises In stability analysis as in steady 
flow ana'ysls: what conditions should be imposed at the separating con- 
tact line? A comparison of the stability predictions for three of the 
postulated boundary conditions is carried out in Chapter 3. The main 
result Is that stability analysis under the lubrication approximation is 
unreliable because the results are sensitive to the choice of boundary 
condition, as .shown in Figure 1.8. 

Finite Element Stability Analysis. The Idea of finite element sta- 
bility analysis Is to subject the finite element estimate of a steady 
flow to all Infinitesimal disturbances of velocity, pressure, and free 
surface that can be represented with the very same set of basis functions 
jsed to estimate the steady base How. Three-dimensional disturbances 
ire represented by finite element basis functions multiplied by Fourier 
tasls functions, i.e. sines and cosines: 

u' » £ u^ 1 (x,y).[i_[ ws(Nz)+jJ cos(Nz)+fck s1n(Nz)]e" U (1.5.1) 
P' ' J Pj^fx.y) cos(Nz)e* u * {1 , 5 . 2) 
h' = J hjp'fx) cos(Nz)e" u 
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Here the uj. pj and hj are nodal parameters that determine the distur- 
bance node; ♦'(x.y), f'u.y) and c '{ x ) are finite element basis func- 
tions; K Is the transverse wavenumber of the disturbance; and the expo- 
nential decay rate x determines shether a . dl sturbance decays, in which 
case X is positive, or grow. In which case X is negative. The analysis 
leads to a large algebraic eigenproblem of the form: 

Vj^Vj. (1.5.4) 

Entries In the Jacobl an matrix, represented by J 1Jt measure the sensi- 
tivity of the residuals of the NavW-Stokes equations and boundary con- 
ditions, as given In C1.2.iWl.2.7). to changes In the nodal velocities 
and pressures, i.e. the finite element solution. The residuals measure 
whatever momentum and mass discrepancies remain In the solution. Entries 
In the mass matrix M y measure the overlap-of basis functions. A solu- 
tion of (1.5.4) 1s an eigenvector x jf which represents a possible mode 
of instability, and an eigenvalue x, which Is the decay rate of that 
•node, i.e. the negative of Its growth rate. In contrast to the hydro- 
static stability analyses carried out by Brown { 197 9 ), i„ wh 1cn ^ elgen- 
problem is symmetric, .It is asymmetric for film-flow problems. This is 
a crucial complication because methods to solve asymmetric eigenproblems 
sre still in their infancy. 

A general method for determining the stability of a coating flow is 
laid out in Chapter 5. The same strategy. Is used to represent pertur- 
>at1ons to a steady two-dimensional base How as illustrated in Figure 
-.3. I.e. the How domain is divided into three zones. Perturbation 
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analyses worked out in Chapter 2 describe the three-dimensional flows (n 
the upstream channel and downstream free surface zones. Matching con- 
ditions at the Inflow and outflow boundaries of the forming zone result 
in Robin conditions similar to the one In Equation (1.4.1). Only the 
most dangerous disturbance modes. I.e. the eigenvectors which satisfy 
(1.5.4) and correspond to the fastest growth rates, are of interest. 

Smart's method Is used here to track the most dangerous distur- 
bance modes. It Is described in Chapter S and more fully In Appendix A. 
Stewart's method Is a comparatively new and sophisticated algorithm based 
on an old sure-fire method known as inverse subspace iteration, which 
relies on repeated solutions of a set of linear equations. Ordinarily, 
complex arithmetic f, required because the eigenvalues and eigenvectors 
of an asymmetric eigenproblem are generally complex*. However, the need 
for complex arithmetic is cleverly avoided by the use of a special basis, 
the quasl-Schur basis, to represent the nost dangerous disturbance modes. 
Because real arithmetic is used exclusively, computational efficiency 1s 
Improved and computer memory requirements are reduced. Furthermore, the 
rate of convergence is enhanced by a Schur-Rayleigh-Ritz scheme, which 
brings about an Improvement in the basis defining thedominant distur- 
bance modes. Continuation. I.e. using a previously calculated set of 
dominant disturbance modes as initial guess. Is used in the implementa- 
tion of Stewart's method to cut computational costs. The overall scheme 
is efficient: so efficient that It usually costs only a little more to 
make one stability calculation than to generate a steady two-dimensional 
flow field. 



2* 

The general method of finite element stability analysis laid cut in 
Chapter 5 1s Illustrated in Chapter 6 with exponential slot and knife 
coating. The main findings are how coating stability depends on each of 
the operating parameters, and ais~ *"iy *h*. Ebbing instability arises. 
The ribbing instability is she- >:-at( *y Increasing capillary 

number Cas uU/c, decreasing coar - . -r increasing the separating 
contact angle that the meniscus ae die, decreasing channel 

divergence angle, or dec^ds ig Reynolds number Re = Ud/v. The stability 
predictions are In close agreement with the results of the visualization 
experiment pictured in. Figure 1.7. The close agreement appears to vali- 
date the accuracy of the Galericln/flnlte element approximation to the 
stability problem. A comparison in Figure 1.9 shows that the lubrica- 
tion stability predictions are off the mark.* Evidence that the cause of 
failure is due to the breakdown of the assumption of nearly rectilinear 
flow In the region near the meniscus 1s presented in Chapter 6. 

The method of finite element stability analysis described 1n Chapter 
5 Is especially useful because 1t not only predicts the stability of a 
flow, it also points out the node of instability. Knowing the mode of 
instability Is essential to understanding the mechanism . I.e. which 
forces contribute to instability and where the seat of instability Is. 
The mechanism of ribbing Instability in slot and knife coating Is exam- 
ined 1n Chapter 6. .The main conclusions are that normal viscous stresses 
acting on the meniscus near the contact line play an important role in 
the mechanism of ribbing, and that the ribbing Instability 1s largely 
explained by examining the normal tractions on the meniscus near the 
contact line. 
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CHAPTER 2 - STEADY. THREE-OIMCHSTQHAl. LI0uid „ lm B E DISTRIBU ti Q h 0H , 
MOV IMG WEB " " 

2.1. Introduction 

Oownstream of almost every coating device Is a zone In which the 
deposited film approaches an asymptotic regime. The asymptotic regime 
is usually solid body translation or nearly solid body motion because the 
force of gravity. Is small In comparison with the viscous drag exerted by 
the coving «b. and furthermore the' role of gravity vanishes altogether 
when the web is horizontal. In this chapter the steady approach to a 
downstream plug-flow regime, I.e. . solid body translation. Is analysed. 
The analysis Is of the decay of small, three-dimensional departures up- 
stream from the downstream regime when the. deposited film Is incompres- 
sible Newtonian liquid having uniform viscosity and surface tension. 
Gravity is regarded as negligible, and thus either capillary force, vis- 
cous force, or. both lead to downstream decay of an upstream departure 
from plug flow. Viscous force resists liquid deformation and thus plays 
a dual role: (I) it opposes transverse. I.e.. cross-web redistribution 
of liquid and (2) it causes decay of filmwlse velocity gradients. Cap- 
lllary force resists departure from planar meniscus shape, and thus It 
*Uo plays a ml-rolt: (I) It promotes transverse redistribution of 
liquid that tends to flatten the meniscus and (2) it counteracts stream- 
wise meniscus curvature which results from decay of filmwlse velocity 
gradients. Three possibilities are Illustrated In Figures 2.1-2.3. 
Figure 2.1 shows a transversely symmetric (In the z-dimenslonj distur- 
bance decaying to a plug profile. Here viscous force is driving and 
capillary force is resisting decay to uniform plug flow. Figura 2.2 
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shows the situation when a coating device 1s either not well aligned or 
nonunifonnly fed: plug flow Is reached only after the liquid fs redis- 
tributed across a long transverse distance. Here the primary effects of 
viscous and capillary forces are. respectively, to resist and to drive 
the transversa redistribution of liquid. In Figure 2.3 the coating Is 
Hbbed by a fluid mechanical Instability In the fornring zone; here, the 
liquid need redistribute only across a relatively short transverse dis- 
tance. In this situation the roles of viscous and capillary forces are 
truly dual because fluid redistributes In the downstream and transverse 
directions over approximately the same length scale. This chapter focu- 
ses on understanding the rate and manner of spatial decay of a disturbed 
flow to Its final plug flow state. There are two primary applications: 
the firs* 1s coating die design and the second Is coating stability. 

While the object of most coating devices Is to apply a uniform layer 
of liquid onto' a moving web. In practice the liquid film is marred by 
Irregularities due to coating die Imperfection or misalignment, mechani- 
cal vibrations, pressure nonuniforaity or fluctuation,, and perhaps even 
fluid mechanical Instabilities - such as the ribbing Instability. What 
<s crucial Is that these Imperfections be held down so that the final 
film thickness 1s maintained uniform to within some tolerance. Surface 
leveling may play an Important role In achieving this goal by reducing 
the -npiitude of coating Imperfections. Quantitative estimates of 
leveling rates could aid significantly In the design of coating dies by 
making clear which disturbances level so slowly that they cannot be 
tolerated. 
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This analysts also applies to the construction of a Robin boundary 
condition at the matching plane between forming and free surface zones 
(see Figure 5.1). The Idea Is that knowledge of how a How decays to 
the final plug regime can be Incorporated advantageously Into a boundary 
condition at the outflow boundary of a forming zone. The advantage of 
using a Robin boundary condition over Dirichlet or Neumann conditions' at 
the outflow boundary downstream of an exponential slot coater 1s demon- 
strated in Chapter 4. Hlggins 1 (1982) analysis is sufficient to con- 
struct a Robin condition for two-dimensional Hows like the one consid- 
ered in Chapter 4 t but a three-dimensional asymptotic analysis' is neces- 
sary whenever three-dimensional Hows or the stability of twc^dlmensional 
flows to three-dimensional perturbations are sought. The asymptotic 
analysis of this chapter is particularly useful for Investigating the 
ribbing instability, the subject of Chapters 5 and 6. 

Several papers have dealt with the onset of three-dimensional rib- 
bing instabilities 1n thin liquid films (see for example Pearson I960, 
Pitts and Grelller 1961, Taylor 1963, and Savage 1976), but to the best 
of the author's knowledge only Fall (1978) has Investigated the down- 
stream decay of ribbing disturbances. Fall used lubrication theory to 
analyze the problem. His results are limited to small wavenumb$r and 
large capillary number (Fall's claim that his analysis Is valid for 
small capillary number Is shown to be incorrect in Section 2.3). 

Much more headway has been made in understanding the approach to a 
final flow regime *hsre upstream departure from the final regime Is two- 
dimensional rather than three-dimensional. Higgins (1980, 1982) was the 
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first to treat the two-dimensional problem In Its full generality. He 
even considered the effect of gravity, but could solve the resulting 
equations only when either gravity or convective Inertial forces could 
be neglected. The mathematical analysis used here 1n Section 2.4 resem- 
bles closely that of Hlgglns. Hlgglns 1 analysis fn turn was Inspired by 
that of Wilson (196*9), who examined the approach to Polseullle flow 1n a 
channel. Both Hlgglns and Wilson were able to use a streamf unction rep- 
resentation for thelr two-dimensional flows; here the three-dimensional 
analogue of streamf unction, vector potential (see Ar1s 1962 and Horse and 
Feshbach 1953) Is used to facilitate representation of three-dimensional 
flows. 

Earlier investigations that wereless general than Hlgglns 1 Include 
those of Bretherton (1961). Cox (1962), Coyne and Elrod (1969), Grcenveld 
and van Dortmund (1970), and Ruschak (1974). Both Bretherton 1 s and Cox's 
analyses were aimed at understanding their experimental Investigations of 
bubbles moving through capillary tubes. Coyne and El rod were Interested 
In the final acceleration of liquid In the cavltated region of a journal 
bearing. Hlgglns (I960) pointed out an error In one of Coyne and Elrod's 
boundary conditions, however, thereby showed their results to be 
Incorrect. Groenveld and van Dortmund examined the liquid film entrained 
on a web being withdrawn from a bath, I.e. dip coating. Nona of these 
analyses account for liquid Inertia, and so represent no more than lim- 
iting cases. Ruschak investigated the same problem as Hlgglns, but his 
analysis Is valid solely for very large surface tension as well as neg- 
ligible inertia. 
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Orchard (1962) and Anshus (1973) Investigated the problem of the 
temporal leveling of small two-dimensional and three-dimensional initial 
disturbances of a stagnant pool, respectively. Both analyses are lim- 
iting cases of negligible Reynolds number, I.e. the acceleration term in 
the momentum equations Is set to zero. 

Two methods are used here to analyse the asymptotic approach to plug 
How. The first begins with the lubrication assumptions and leads to an 
approximate solution, valid -hen inertia Is small and disturbances have 
long transverse wavelength. I.e. small transverse wavenumber. The second 
method employs no assumptions other than that disturbances be small, and 
leads to a result which Is valid for the entire range of parameters. The 
second method Is to construct a solution of the exact governing equations 
linearized about the final plug flow regime. The linearized equations 
are of Oseen type (see Oseen 1927). 

Both the lubrication and Oseen solutions are constructed with the 
aid of normal mode analysis. Arbitrary disturbances can be represented 
as Fourier series or Integrals - depending on whether tie finlteness 
of the web In the transverse direction Is taken Into account or not. 
Because the lubrication and Oseen equations are linear, disturbances of 
Individual wavenumber can be singled out and studied Independently, a 
solution could then be constructed for a given upstream disturbance as 
a series, or Integral In wavenumber. However, the construction of such a 
series or Integral solution would be difficult in practice. Fortunately, 
much can be learned by studying the decay of disturbances which have 
single Fourier components. Furthermore, linearized three-dimensional 
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stability analysis ts performed by Investigating one Fourier component 
at a time, and so it Is desirable to do the same In constructing a Robin 
boundary condition (see Section 2.6 and Chapter 5). 

Both the lubrication and Gseen analyses reduce to eigenproblens. 
Each solution Is an elgenpalr, i.e. an eigenvalue and corresponding 
e1genfunct1on k and represents a possible node of disturbance and Its 
spatial decay, rate. All disturbances decay with downstream distance 
because no forces which could cause disturbances to 3 row with downstream 
distance are Included in this analysis. The characteristic equation to 
which the Oseen analysis leads has an Infinite mmber of roots. Any 
disturbance of specified wave number could be represented as a series In 
the corresponding elgenfunctlons. In general, some 0 f the eigenvalues, 
I.e. spatial decay rates, are complex. These correspond to modes which 
are oscillatory in downstream coordinate. J.e. damped standing waves. 
Real eigenvalues on the other hand correspond to cades that decay expo- 
nentially In downstream distance. By contrast with the Oseen problem, 
the characteristic equation of the lubrication analysis has only two 
allowable roots, and both are real. Physically, this is because the 
luortcation assumptions dictate that velocif disturbances be parabolic 
in the fUmwIse coordinate, y. Thus, the possible disturbances are 
severely limited. 

In Section 2.2 the Oseen equations, boundary conditions, and the 
lubrication equation are laid out. In Sections 2.3 and 2.4 solutions 
of the lubrication and Oseen equations are constructed, respectively. 
Section 2.S contains limiting and numerical results for the Oseen solu- 
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tlon. The Oseen and lubrication results are discussed and compared with 
those of Fall and Higgins in Section 2.6. Finally, in Section 2.7, 
Wilson's (1969) analysis Is extended to include three-dimensional dis- 
turbances. 

2.2. Oseen and Lubrication Equations 

The governing equations and appropriate boundary conditions for vis- 
cous free-surface flows of incompressible, isothermal. Newtonian liquids 
were given in Equations (1.2.1 J-(1.2.3) and U.2.5J-U.2.7). In the 
free surface zone downstream of a_fonr.ing zone, the linearized form of 
these equations applies. The linearization is about the final uniform 
plug flow state: 

' «•!. P'-O. N-l , (2.2.1) 

where U is the velocity. P is the pressure, and H is the elevation of 
the free surface above the moving web. The capital letters designate 
the base flow. The governing equations linearized about this state are 

M Reif a ' 2 « -*P • (2,2.2) 

7, i 8 ° • . (2.2.3) 

where s UH/v, u being the web speed, v the kinematic viscosity, u 
the velocity, with components (u.v.w), and p the pressure. The. lower 
case letters designate a disturbance to the base How. The disturbance 
must obey the following set of linearized boundary conditions: 
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^ 8 f - 0 . (2.2.4) 

(t-i)xi-0 • y-1., (2.2.S) 

N cai : il" 7 ?l h »yl. (2.2.6) 

A'l'S ■"• y-1 . (2.2.7) 

Equations (2.2.4)-(2.2.7) are the no-slip and- no-penetration, shear 
stress, normal stress, and kinematic conditions, respectively, for a 
small disturbance from unlfono plug flow, t is the stress tensor, which 
under the restrictions that the liquid Is Incompressible and Newtonian 
Is given by 



t = (7u ♦ (7u) T ] - 7p , 



(2.2.6) 



M ca s uil/a Is the capillary number, u being, the dynamic viscosity and o 
the surface tensior., and is the two-dimensional gradient In the 
plane of the moving web. so that i S the two-dimensional Laplaclan. 
Disturbances must eventually decay to zero at very large distances down- 
stream; hence the additional restriction that 



u ♦ 0 , p v o , and h ♦ 0 as x ♦ 



(2.2.9) 



must be satisfied by any allowable disturbance. The disturbance velocity 
must satisfy the matching condition 

JLhjM) -u^y.z) 6 x = 0 . (2 . 2 . 10) 
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where ^(y.z) i s the velocity disturbance at the upstream boundary of 
the free surface zone (see Figure i.i). The solution of Equations 
(2.2.2)-(2.2.lO) is constructed In Section 2.4. 

Equations (2.2.2M2.2.10) can be simplified by means of the lubri- 
cation approximation (see Cameron 1976, Mitchell 1950, andTlpel 1962), 
which rests on the following conditions: (I) convoctive Inertlal forces 
are negligible and thus Reynolds number, M Re , Is effectively zero; (2) 
flow Is nearly rectilinear, and hence the fllmwlse velocity component, 
v. 1s small compared wlth.streamwlse and transverse velocities, u and 
w. respectively, and fllwlse gradients are larger than streamwlse and 
transverse gradients; and (3) gravity Is negligible (already Imposed 
above), with these assumptions, Equations (2.2.2)-(2.2.10) are reduced 
by the standard scaling arguments of lubrication theory to: 

3 2 u„/Jy 2 « 7 n p . 

7 • U ■ 0 

3u n /3y = 0 

J_ • u 3 0 
£ u * 0 , h ♦ 0 
u n {x,y,z) * ^{y.z) 





(2.2.11) 




(2.2.12) 


0 y a 0 


(2.2.13) 


9 y = 1 


(2.2.14) 


9 y » I 


(2.2.15) 


9 y » 1 


(2.2.16) 


as x • • 


(2i2.17) 


9 x » 0 


(2.2.18) 
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Here the subscript II again Indicates that the vectors have components 

only the streamwlse and transverse directions, x and z. respectively. 
The pressure field can be found from equations (2.2.11)-(2. 2 .i4) In the 
usual way (see for example, Catron 1976) by solving Equation (2.2.H) 
with boundary conditions (2.2.13) and (2.2.14), integrating that result 
and Equation (2.2.12) over y fro* the web to the free surface with 
Equations (2.2.13) and (2.2,16) „ boundary conditions, and combining 
the results of. these integrations to eliminate velocity. The result is 
'il ' (h 3 v n p/6) - Oh/3x)/2.. ( , 219) 

This differs from the well-known lubrication result by a factor of one- 
half on the left side because the velocity profile is semi -parabolic 
rather than parabolic. Finally, because pressure I, not a faction of 
the fn«lse coordinate, y. Equation (2.2.^5) can be used to eliminate 
pressure In favor of film thickness: 

7 II ' ("Si^hlJ = 3« M aW? (2 2 2Q) 

The gradients of free surface elevation, h. m small by assumption, and 
so It is consistent to linearize about the base state. H = 1. Therefore. 
Equation (2.2.20) Is approximated by 
4 

fj.h = 3N ah/ax . 
11 ca • • (2.2.21) 

To the author's knowledge, this equation has not appeared in the litera- 
ture before. Solutions are constructed in the following section. 
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2.3. LubHcatlon Solution 

Solutions of Equation (2.2.21) can be found by separation of vari- 
ables and have the form 

h ■ t(x) cos(Nz) + ?(x) sin(Nz) , (Z.2A) 

where M is the wavenumber and Is inversely proportional to the wavelength, 
which is 2t/M. Inserting this in Equation (2.2.21) produces the ordinary 
differential equation: 

[t"" - 2N 2 r« ♦ 3N ca r' + N 4 t]cos(Nz) 

♦•[?- - a! 2 ?" * 3M ca ? ♦ N 4 ?]sin(Nz) -0. (2.3.2) 

Both bracketed terms oust be equal to zero because z is arbitrary, and 
thus 

r- - 2N 2 t" * 3M ca T' > H\ a 0 5 t»"- 2H 2 ?" ♦ 3M f * N 4 ? » 0 . 

Ca 

(2.3.3) 

Because both t and ? satisfy the same differential equation, It suffices 
to solve for t alone. The auxiliary equation for (2.3.3) is 

where t * ce' ax and a Is the spatial decay rate. Only roots a of 
Equation (2.3.4) which have positive real part are of interest because 
otherwise the requirement that disturbances disappear far downstream, 
Equation (2.2.17), would be violated. Descartes' rule of signs shows 
that Equation (2.3.4) has either two or no positive real roots. By 
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rewriting It. in the form 

((«/N) 2 - l] 2 - 3 (Nca/N 3 ) (o/M) » 0 , (2 . 3 . s , 

1t can be seen that («/N) depends solely on (N ca /N 3 ), i.e. the ratio of 
exponential decay rate to wavenumber is a function exclusively of the 
ratio of capillary number to wavenumber cubed. Equation (2.3.S) can be 
solved by perturbation analysis (see Lin and Segal 1974) when the param- 
eter (N ca /N 3 ) l s small; the result Is 

In the other limit, that Is when (N ca /N 3 ) Is large, perturbation analy- 
sis leads to the roots 

a/N " > (3H ca /H 3 ) 1/3 . (2 . 3 . 7) 

The first result In Equation (2.3.7) matches Fall's (1978). The solu- 
tions at intermediate values of the parameter (N CJ /N 3 ) are easily found 
by Newton's method. The smallest positive root, (a/N). of Equation 
(2.3.S) for parameter (N ca /N 3 ) ranging from C.01 to 10.0 is shown in 
Figure 2.4. 

2.4. Oseen Solution 

Any solenoidal vector field, i.e. one that satisfies the continuity 
equation (2.2.3), can be represented by a vector potential (see. e.g. 
Arls 1962): 

JL*'xA . V.A-0. (2<4>1) 
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The divergence of the vector potential A does not affect toe velocity 
field, u, and therefore can be set to zero as In Equation (2.4.1). This 
property, that the divergence of the vector potential can be chosen 
arbitrarily, 1s known as gauge imparlance (Morse and Feshbach 1953). 
Taking the vector potential to be sol'enoidal simplifies the curl of the 
velocity field: 

' x u » V x V x A « - v 2 A ♦ V(r .£) » - 7 2 A (2.4.2) 

which Is needed below. The pressure field can be eliminated from Oseen's 
equation. Equation 2.2.2, by taking the curl to get an equation for vor- 
ticlty, i.e. the curl of the velocity field. Equation (2.4.4) can then 
be used to replace vortldty by the vector potential to give 

- "Re |7 ' 2 i ■ £ • (2.4.3) 

This equation has a separation of variables solution of the fora 

A«e" ax [a(y) sin(Nz) +i(y) cos(Nz)] , (2.4.4) 

where H fs again the wave number in transverse coordinate, z, and * Is 
exponential decay .rate in the streamwise coordinate, x. Replacing A In 
Equation (2.4.3) by Its equivalent In Equation (2.4.4) and dividing out 
the consnn e" oX gives the fourth-order differential equation: 

(a- ♦ [2(aV) + aH^Ja- ♦ [(a 2 -N 2 ) (, 2 -M 2 + ,M Re )]aj sin(Nz) 
* {1"" ♦■ (2<«V) ♦ oN Re ]a-> [(« 2 -M 2 ) (« 2 -N 2 ♦ aN Re )Ji}cos(Nz) = 0, 

(2.4. S) 
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a prime denoting differentiation with respect to y. Because the origin 
of r Is arbitrary, each bracketed expression must be zero: 

i" ♦ [2(«V) ♦ , NRe Ja- ♦ [(.V, {8 2. M 2 , ^.o, (2 . 4<6) 

and likewise fori. Because the differential equation for i is Iden- 
tical, it has the same general solution as (2.4.6), namely 

a * Cj cos(ay) + Cj sin(ey) + c^cosfry) + c, sln(ry) = c f S f 

(2.4.7) 

where 8 = « - N 2 . Y 2 5 a 2 - n 2 + «N Re> £f Is a vector containing three 
arbitrary constants, and 

[S,J 5 [cos(sy), sin(sy>, costry). sin( Y y)] ( 2 .4.8) 

is matrix vector that Is convenient for subsequent manipulations. Sum- 
mation over a repeated index is issued. Each of the vector coefficients 
c f contains three scalar coefficients, for a total of twelve, a also 
contains twelve scalar coefficients, and thus by (2.4.4) the vector po- 
tential contains twenty-four. However, the requirement that A be sole- 
noidal in (2.4.1) imposes eight linear constraints on the twenty-four 
unknown scalar coefficients. Furthermore, four of the scalar. coeffi- 
cients do not affect velocity field. To understand why this is so. it 
Is instructive to express the vector potential in Helmholtz' representa- 
tion (see e.g. Arls 1962): 

-here * is a scalar function and * is i'a solenoidal vector field. The 
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requirement that Abe solenoldal implies that* must satisfy Laplace's 
equation: 

V . A » 7 2 $ = 0 

There are four solutions of this equation that have the same form as 
terms 1n A, which are 

♦ 3 ti x cosey + * 2 s1nsy)e" ax cos(Nz) • 
♦ {* 3 cosey +.* 4 s1nsy)e" aX s1n(Nz) . 

This result 1s not fortuitous because the differential operator In Equa- 
tion (2.4.3) for the vector potential can be factored into two commuting 
operators, and thus A * Aj + where 

^'0 , (V 2 • i-jAg » o . 

The first of these Is the vector Laplace equation. . It is therefore 
possible to find solutions * which both satisfy Laplace's equation and 
can be used 1n Helmholtz' representation of the vector potential. But 
their contributions to the vector potential disappear from the curl of 
the vector potential, I.e. the velocity field. 

The final cc»nt of Independent scalar coefficients in the vector 
potential A that affect the velocity field is twelve: eight of the ori- 
ginal twenty-four scalar coefficients are restricted by the requirement 
that A be solenoldal (2.4.1) and four more play no role in determining 
the velocity field. To find the remaining twelve unknown coefficients 
would require that boundary conditions (2.2.4)-(2.2.7) be translated 
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into ones on a and a. This Is difficult because pressure cannot easily 
be eliminated from the normal stress equation (2.2.6). it -is much more 
convenient to return to the primitive variables, u, p, and h." u can be 
constructed by talcing the curl of the vector potential A and has the 
following form: 



u 2 u 3 



u + u - 



|"icos(N2)-j |- Ul 
I i cos(Nz) I I vj v 2 v 3 



•1 



Ljc cos(Nz) J w 2 * 3 w 4 j 



p_ cos(Nz)l 
I J.cos(Nz) 
l_lc cos(Nr)_ 



r 5 r 5 2 



v 2 v 3 v 4 



w 3 *4- 



pos(sy)*j 

I stn(sy) I 
I i 
I cos(yy) | 

Lsfn(Yy) J 



-OX 



(2.4.9) 



Here 0 t 5 x cos(Nz), Og 5 £ cos(Nz). Og s k sin(N 2 ), fl ; s L $ 1n(Nz), 
5g s i sln(Mz), and ^ » k cos(Nz) are vectors and ufj and are the 
elements In the 1-th row and J-th column of the velocity coefficient 
matrices which appear in Equation (2.4.9a). Coefficients u{j and ufj 
are related to c, and c f by the general solution for veloc1ty.(2.4.9), 
the relationship between velocity and vector potential (2,4.1)', and the 
equations that give the general solution of the vector potential - 
(2.4.4). (2.4.7). and the equivalent of (2.4.7) for a. This relation- 
ship is displayed In Figure 2.5. where and ii fj are transformation 
matrices that connect velocity coefficients with vector potential coef- 
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flelents. Each transformation ntrtx has t»., null vectors as shown, and 
therefore only twenty of the twenty-four velocity coefficients are Inde- 
pendent. The null vectors result from contributions to the vector 
potential which are gradients of a scalar function, as explained in the 
preceding paragraph. Only eight of the twelve scalar coefficients in 
the four c, are Independent, as explained above, and likewise for c f . 
Thus, only six of the twelve scalar coefficients ufj are independent, 
and likewise for tf y Pressure can be found by substituting velocity 
as given by Equation (2.4.9) into the Oseen equation (2.2.2) and inte- 
grating. It has the following form: 

posfflyrj 

P ♦ P ■ (cos(Nz)[ Pl .p 2 ,p 3 ,p 4 j ♦ s1n(Mz)[p l .p 2 ,p 3 ,^J } , IStn(By) ie-« 

Icos(Ty)l 

JjlnfyylJ 

h (cos(Nz) Pl ♦ sinfMzJp^S, , (2 . 4>10) 

where p, and p | . the pressure coefficients that appear in (2.4.10a). are 
related through the Oseen equation to the velocity coefficients in 
(2.4.9). and thus are not Independent. The free surface elevation can 
be shown by substituting velocity from (2.4.9) Into the kinematic bound- 
ary condition (2.2.7) and Integrating In the streamwise coordinate, to 
have the form: 

h * "h » (cosOlzlh! ♦ slntNzJh^e-" , ' [ZAAl) 

where h L and hj are scalar unknowns.. 
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The simplest means to find the linaar constraints on the velocity 
and pressure coefficients defined in Equations (2.4.9) and (2.4.10). is 
to return to the basic equations that they must satisfy, namely the 
Oseen and continuity equations, (2.2.2) and (2.2.3). Before doing this. 
It is convenient to choose a particular phase of the perturbation in 
transverse coordinate, z: a crest of the free surface is located at 
z • 0, i.e. 3h/3z » 0 and aV»z 2 .< 0 there, With this choice, the 
coefficients that have tildes in (2.4.9)-(2.4.ll) drop out. (Were the 
choice to retain only the coefficients with tildes, the following 
results would be Identical except' that the signs of some terms con- 
taining wavenumber would change; the characteristic equation for decay 
rate « would be Identical.) The linear constraints on ^ and Pf from 
the Oseen and continuity equations are 

{ *ti y \ * R Jk) * "PA a 0 - 

(u 2^Sk * R jJ -PrjlA" 0 (2.4.12) 

( - M ?j * U 21 R 1J + «"!jJSj • 0 (2 .4. 13) 

where 6 Jk is the Kronecker delta, i.e. « jk » 1 1f j = k and 5jfc = 0 
otherwise, and 
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r° " fl 0 °i p 2 o o o -i 

^.|:.::;i.K J)i(w --i:: 2 '°:i 

L° 0 Y °i lo 0 0 Y 2 



Y J 
(2.4.14) 



Because each S, is a „ independent function of the mnwfse coordinate y 
(see Equation (2.4.8) ) and Equations (2.4.12) and (2.4.13) mst hold 
throughout the domain, it follows that: 

U 3jt Y Sk + R J k J + »P|c °° . 

* ° U U * U 2 C 1 R 1J + " 0 J- 1,2,3,4. (2.4.16) 

These are sixteen equations In the sixteen scalar coefficients of velo- 
dty and pressure. However, only ten of the equations are ^dependent, 
leaving six coefficients unknown as anticipated on page 2-17. Ten inde- 
pendent equations from the set of sixteen are: 

a *«e u l + -Pi. - 0 ' (2 . 4 . 17) 

a " R e "2 * "P 2 = ° . (2 . 4 . 18) 

Re 1 P 2 (2.4.19) 
Re 2 P l (2.4.20) 



p 3 = 0 
P 4 3 o 
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aN *e W * +Hp l 3 ° (2.4.811 

"* W 2 * NP 2 ' ° (2.4.22, 

(2.4.23) 
(2.4.24) 

■--j'YV^-b ( 2 .4. 25 , 

The Oseen solution 1s completed by enforcing the boundary conditions 
at the moving web. at the free surface, and at Inflow and far-downstream 
boundaries. The no-sllp and no-penetrat1on boundary condition (2.2.4) 
Imposes three constraints on the velocity coefficients: 

Ui ♦ u- » 0 

1 • 3 (2.4.27) 
v, *■ v, « 0 

13 (2.4.28) 

W, ♦ W, »' 0 

13 (2.4.29) 

The kinematic condition (2.2.7) relates tie free surface elevation per- 
turbatlon to velocity: 

+ ,2.4.30) 

The shear stress boundary condition (2.2.5) dictates that both components 
of viscous shear at the free surface must be zero: 

- « fed) ♦ u| iVj (l> - 0 {2<4>31J 
»{ f Vj. ( » - oU 21 S i'H =0 (2.4.32) 
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The normal stress boundary condition (2.2.6] equates, the viscous normal 
stress and pressure at the free surface with capillary stress due to., 
meniscus curvature: 



2u 2 V 1J S J (I) - Pj S t (l) .J_ 6 2„ ia0 . (2.4.33, 

ca 



Equations (2.4.17)-(2.4.33) are seventeen homogeneous scalar equations 
for the seventeen unknown coefficients - twelve in the velocity field, 
four in the pressure field, and one fn free surface elevation, as defined 
1n Equations (2.4.9M2.4.11) - plus the unknown decay rate. a. Equa- 
tions (2.4.17)-(2.4.33) appear In matrix form In Figure 2.6, and the 
determinant of that matrix must be zero In order for a solution to exist. 
The condition that this be so Is the characteristic equation for spatial 
decay rate, a. The saventeen scalar coefficients in the velocity, pres- 
sure, and free energy surface elevation can be determined .once a Is 
known, but only up to an arbitrary constant. 

It 1s convenient to eliminate most of the unknown coefficients in 
the set of seventeen In order to reduce the order of the determinant and 
thereby simplify calculation of the characteristic equation. Equations 
(2.4.17W2.4.29) each contain at most three of the unknown coefficients, 
and so they are easily rearranged and combined to express each of the 
unknown coefficients as functions of a set of four coefficients. The 
optimal set of coefficients found by applying the strategy of Stadtherr 
et al. (1974) are Uj. Ug. and u 4 , and the resulting expressl ons' are 

p l a ' \s "l (2.4.341 
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p 2"-"Re«2 


(2.4.35) 


P 3 a 0 


(2.4.36) 


P 4 " o 


(2.4.37) 


u 3 " - u l 


(2.4.38) 


3 - Bu 2 /a 


(2,4.39) 


v 2 ■ SUj/o 


(2.4.40) 


¥3 3 SU 2 /o 


(2.4.41) 


V 4 3 - 8 2 Ul /aY 


(2.4.42) 


3 Niij/a 


(2.4.43) 


W 2 3 Nt V° 


(2.4.44) 


w 3 3 - Nyo 


(2.4.45) 


W 4 3 (au^ + BYU 2 /a)/N 


(2.4.46) 



All but four coefficients - Ul , u 2 , v and hj - can now be eliminated 
from Equations (2.4.30)-(2.4.33): 

ursine - £ siny] ♦ u^cosy - cosb] ♦ h^sl] = 0 (2.4.47) 

2 2' • 
"lt^y- slnr - £ sine] ♦ u 2 [£ (2cosb - cosy)] + uJcosy) 3 0 

2 2 '(2.4.48) 
slnr - 2 i sinfl] * u 2 (i (cosb * JL±- cosy)] ' 

Y 2 Y N Z 

+ u 4 [-rCOSY] « 0 (2.4.49) 
M 

ufiicou - cosy) ♦ ^21 cos6] + ^[ 2 ( IJflB ' . I sinr) + f!!|ej 

+ h ltr"i * 0 8 (2.4.S0) 
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Equations (2.4,47) and (2.4.501 can be combined to eliminate h L and 
Equations (2.4.48) and (2.4.49) to eliminate u 4 . The characteristic 
equation Is found by setting the determinant of the two resulting equa- 
tions equal to zero: 



2 2' 

- -~- slny ♦ 2s1nfl 
Sy 


1 
1 
1 


- 2cos6 ♦ 


8 


S 2 ft 

N ca Y 


1 
1 


jji{COSY 

N ca 


- C0S8) 


2 2 

M^-^L-COSB - 2s cosy) 

D 


1 
1 


' 8 


slnfl - 2y slny) 



* 0 
(2-4.51) 



In the limiting case that wavenumber equals zero, this equation Is equiv- 
alent to Higglns' (1980) Equation (2.4.40). When wavenumber, N, equals 
zero the following simplifications result: (I) 6 is equal to a, the ex- 
ponential decay rate and (2) Y 1s equal to /a* + a N Re . These quantities 
correspond to Higglns' 1n the following w*y: (1) a used here Is the 
negative of Higglns' a n (compare Equation (2.4.4) with Higglns 1 Equation 
(2.2.21)); (2) therefore, 6 Is also the sane as the negative of Higglns 1 
<* n ; and (3) y Is the same as the negative of Higglns 1 x^. When Higglns 1 
symbols are subsltuted for a, b, and y, the first column in Equation 
(2.4.51) is scaled by y/B, and the two columns are Interchanged, the 
result is Higglns 1 Equation (2.4,10). Equations (2.4.51) can be 
rearranged and simplified to get the result: 
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^8 [s1n( 8 )cos(Y) -ieotdhinhj] ♦ (5 S WyV)cos(s)cos( y > 
♦ i (6 4 ♦ 66 2 r 2 * T*j3ln( B )stn(T) - 4sV ♦ y 2 ) * 0 \ (2.4.52) 
This equation is transcendental in decay rate, .. The coefficients, „, 
s«d can be determined up to an arbitrary constant from the elements' 
of a col.,, of the adjugate (see toundson 1966, who uses the word adjoint 
In place of adjugate) of the matrix in Equation (2.4.51): 

"l-«»»-^'«(r>' ,(2.4.53) 
2 2 

u 2 = 2s1n(B) -I_tS_ s1n ( Y) 



~ 5iniT ' ' (2.4L54) 



The coefficients, ^ and u 4 are found. by Inserting the above equations 
for u : and ^ in Equation (2.4.47) and (2.4.48), respectively: 



h l a -?ffsin{ 8 )cos(T) -icos(a)s1n( Y )J .(2.4.55) 



u 4 --2i s i„( B)t I?l|! sfn(y) . {2As6) 

6 

The remaining coefficients are easily found from Equations (Z.^)- 
(2.4.46). 

A unique solution can be constructed by invoking the Inflow and far- 
downstream boundary conditions. The far-downstream boundary condition 
(2.2.9) restricts the allowable values of the decay rate a to toe set of 
c^lex numbers with positive real parts, because those ar« the ones that 
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correspond to disturbances which decay with downstream distance. The 
decay rates which belong to the set of complex numbers with negative, 
real parts correspond to modes which decay with upstream distance, and 
thus represent solutions which have an upstream asymptotic plug-flow 
regime. Examples are uniform film flows which meet stationary obstruc- 
tions, such as the upstream edge of a knife coater. The number of decay 
rates that have positive real part arid satisfy the characteristic equa- 
tion (2.4.52) ,1s Infinite, as demonstrated In Section 2.S. Each decay 
rate corresponds to a distinct mode. I.e. a distinct How pattern. An 
arbitrary perturbation of given transverse wavenumber H can contain any 
or all of these modes. Furthermore, an arbitrary perturbation can have 
Infinitely many Fourier components, I.e. Infinitely many values of trans- 
verse wavenumber H. Therefore, an arbitrary perturbation with arbitrary 
phase 1s described by 

oo <D 

where b tJ and b^ are arbitrary constants, ^ are the wavenumbers which 
fit the transverse dimension, L: 

H i " p ; i s (2.4.58) 

and 0j are the roots of the characteristic equation (2.4.52) that have 
positive real part. Finding the constants in Equation (2.4.57) so that 
u satisfies the inflow boundary condition (2.2.10) would be a difficult 
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(2.4.59) 



task: it would require the solution of the adjoint problem (see Friedman 
1956), u* and u*. so that the biorthogonallty conditions. 

^.aj). u*(H k ,« t )> = V j4 <u(N { , aj ), <L*(N 1 ,a J )> 

^.aj). u*(N k ,« t )> » 0 
^.aj). U*(N k .a t )> » 0 . . 

could be used. Here, the inner product is 
LI 

5 / / (T'b)dydz . 

oo~- 

where the overbar indicates the complex conjugate. If u* and u* were 
known, the constants could be determined by the inner products" 

b lj " < Hotf**»|.«j»>/<u« 1 .« J ). u*(ll i ^ J )> 

*ij 3 %.i*(M 1 .« J W<u(N 1 .a J ). uty-,.^ . (2.4.59) 

However, the difficulties of calculating the adjoint solutions and then 
using then to construct the constants In (2.4.59) can be molded when 
the inflow boundary of the asymptotic domain is far enough" downstream so 
that only one mode, the most slowly decaying one. persists. This simpli- 
fication is used to advantage in the three-dimensional stability analy- 
sis described in Chapter 6 as explained in Section 2.6. 



56 

2.5. Limiting Cases and Numerical Results 

The characteristic equation (2.4.S2) is a complicated transcendental 
function of the decay rate a. Two methods are employed here to solve It. 
The first Is to construct approximate solutions valid only for Minting 
parameter values - the parameters being Reynolds number. N^. capillary 
number. « ca> and wavenumber, N. The second 1s to employ secant itera- 
tion on a digital computer to track the roots of (2.4^52) as parameters 
are varied. The numerical results are set out below in Tables 2.1-2.5. 

A useful limit Is when inertia is negligible compared with viscous 
forces, I.e. when Reynolds number Is approximately zero. However, 
setting Reynolds number to ze;«o 1n the characteristic equation (2.4.52) 
produces the trivial equality 0 s 0. One remedy Is.. to expand each arm 
in a Taylor series In Reynolds number, collect the lowest order terms - 
In this case these are the linear tenns in Reynolds number - and set 
them to zero. The result Is 

2a[cos 2 (8)-6 2 ] ♦ (6 COS(B)s1n(B)-6 2 ]/N - 0 M 0 = 0 . (2.5.1) 
When wavenumber 1s zero this is identical to Higglns' (1980) result 
(compare by replacing both a and 8 in (2.S.1) by the negative of Hlggins' 
;„). Equation (2.5-1) Is not readily solved, but two limiting cases in 
capillary number are easily Investigated by standard perturbation methods 
described in Lin and Segel (1974). 

First, when capillary number is small so that capillary forces dom- 
inate viscous forces the slowest decay rate is 

. 'l**'™^ N c. <<l '«Re 3 ° (2.5.2) 



57 



Thus the slowest decay rate is slightly less than the wavenumber. Thfs 



Is so only when t< ca /M 3 « 1; otherwise 

a l 3 " 4/3N ca ■«"t i «l.Sa'"0. (2.5.31 

so the slowest decay rate goes to zero much more rapidly than wavenumber. 
the results in Equations. (2. 5. 2) and (2.5.3) are Identical to those of 
lubrication theory as given In Equations (2.4.6) and (2.4.7), respec- 
tively. Mote, however, that (2.5.3) is invalid in the two-dimensional 
limit of wavenuniber going to zero because the primary velocity component 
of tlie corresponding mode Is transverse. The easiest way to see that 
this Is so Is to examine the pressure field as given by Equation (2,4.10); 
the streamwlse pressure gradient is proportional to decay rate and there- 
fore goes as the fourth power of wavenumber by (2.S.3), but the trans- 
verse pressure gradient goes as the first power of wavenumber and there- 
fore dominates In the Unit of wavenumber approaching zero.' Because 
pressure gradients are the only forces that drive How when Reynolds 
number is zero, the above assertion that the primary velocity component 
is transverse Is validated. Hlgglns (1980) shows the slowest decay rate 
at zero Reynolds number and small capillary number In the two-dimensional 
limit to be 

«2 3 ' :5M ca il/3 "ca^^N'NRe 30 - (2.S.4) 

More rapidly decaying nodes have decay rats: *h1ch satisfy 

23^ = stn(28, tf2 ) n « 1.2.3,... . (2<5>S) 
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Higglns (198Q) gives the approximate result: 

B 2rr>l 3 (n * t {in[(4«flJ»J ■ = N Rfi = 0 (2.S.6) 

3 2rt+2 3 in +V )] " 4 Zn[(4f * 1,f] "ca = % ' 0 ■ (2.S.7) 



where f s /• ,1 , and a = + u* . 

The second Uniting case at negligible Reynolds number is for very 
large capillary number, f.e. capillary forces small compared with viscous 
forces. The slowest dec^y rate Is 



„ 1 H cosh N s1nh H - H Z 
H " 2 2 

ca cosfT N ♦ n 



ai * """" y;: y " : M ca » 1 * v 3 0 • «-s.8>- 



As before, * x » 0 when « = o, and should be disregarded. Other roots 
satisfy the equation 

cos2fi n =6 n (2.5.9) 
and are approximately (see Higgins 1980) 

6 2 30 ' 73909 (2.5.10) 

6 2wl =■ IB ♦ 1 1n(2r«r) l/H ca = M Re » 0 (2 . 5 . n) 

8 2n+2 - rm - i ln(2n.) W Ci = % * * . (2.5.12) 
where n = 1,2,... . 
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Another useful Halt Is when Inertia dominates over viscous forces, 
i.e. Keynolds' number Is very large. Equation (2.4.52) reduces to 

sln(B)coslr) - icos(8)sln(r) « 0 i/H Re * o . (2.5.13) 

For large wavenumber, I.e. short wavelengths. (2.5.13) can be simplified 
further to get: 

tafl(YlA ' l/M (2.5.14, 
which has roots given approximately by 

%^(nM + H 2 ) M Re »l. (2.S.X5) 

where n » 1,2,... . when wavenumber Is small, the dominant term In 
Equation 12A.SZ) is (SB 2 ♦ 2e 2 r 2 ♦ Y 4 )cos(e)«s(Y>. Setting this to 
zero gives the roots 

l ,Z 9 

Vl "•^ [r(2n * l> 1 m ' *** >>l ' ' • (2 ' S - 161 

Where " * 1,2 dominant decay rate, « l§ Is given by Equation 

(2.5.3) provided N\ e /H ca << 1, which is ordinarily true when wave- 
number, M, is small.. 

The characteristic equation (2.4.52) was solved numerically by the 
IMSL subroutine ZANLYT, which uses secant Iteration (see e.g. Hlldebrand 
1974) to find complex roots of equations. The five slowest decay rates 
were tracked through the three-dimensional parameter-space represented 
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by M Re' M ca' and N by us<n 9 continuation, I.e. by using as Initial guess 
a solution to (2.4.52) at a nearby point in the parameter-space. The 
asymptotic results In Equations (2.5.2)-(2.5.4), (2.5.6)-(2.5.8), 
(2.5.10M2.S.13), (2.4.15), and (2.4.16) provided initial guesses to 
begin continuation and checks to certify that the numerically determined 
roots were indeed the five slowest decay rates. The two slowest decay 
rates, Oj and a n , are given in Tables 2.1-2.5 for Reynolds number 
ranging from 0.01 to 1000, capillary number ranging froa 0.01 to 100, 
and wavenumber ranging from 0.01 to 10, which covers the entire spectrum 
of operating conditions typical of -Industrial and scientific coating 
devices. 

2.6. Discussion 

Much can be learned about the manner and rate of spatial decay of a 
small, arbitrary, steady upstream disturbance In a thin film on a trans- 
lating web from the foregoing analysis. What are of greatest importance 
are the modes which decay the most slowly, because they persist furthest 
downstream. Hlgglns (1982) has described the results when wavenumber Is 
zero. In that limit, Hlgglns (1982) shows that the dominant eigenvalue, 
«j, 1s always real and always belongs to. the sans family of modes, i.e. 
the first and second eigenvalues, a l and <» 2 , never swap positions because 
the Inequality ^ < Real(a 2 ) holds for all values of Reynolds number and 
capillary number. The situation can be much more complicated when wave- 
number is finite. To illustrate the kinds of behavior that can occur, 
Figures 2.7, 2.9, and 2.13 show how the first five decay rates, ^ through 
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<J 5 . vary with wavenumber when Reynolds number Is 100 and capillary number 
1s 100, 1, and 0.01, respectively. 

Manner of Spa tial Decay. When Reynolds and capillary numbers are 
both 100 (Figure 2.7), the first five decay rates - numbered to put the 
real parts in ascending order when wavenumber Is zero - are real and 
retain the same ordering as wavenumber Increases from 0.1 to 10. As 
wavenumber decreases, each of the spatial decay rates except a x approach 
a finite limit. ^ approaches zero for the same reason as In the zero 
Reynolds number case described 1n Section 2.5:- the mode corresponding 
to oj is primarily a transverse redistribution of liquid and as wave- 
number goes to zero the transverse distance over which the film levels 
increases without bound. The higher modes become two-dimensional, I.e. 
they have no transverse velocity component or transverse gradients, In 
the zero wavenumber limit and therefore have finite decay rates. When 
the wavenumber beccses large the decay rates become, large as well because 
the length scale over which liquid must redistribute becomes small. The 
minima In the curves for a 3 and a g are presumably due to coupling between 
transverse and streamwlse velocity components. The real values of spa- 
tial decay rate In Figure 2.7 Indicate purely exponential decay. This 
is Illustrated further In Figure 2.8, which is a portrait of six streak- 
lines of a flow disturbed by the most slowly decaying mode when wave- 
number Is unity/ The streaklines, which are the loci of particles that 
have passed through a given point, were constructed for this steady flow 
by Eulerian time integration of velocity: 



0.1 t ~ 

Wavemooer, N 

Ptqure 2.7. Behavior of the first Ave spatial decay rates 
as uavsnuaber varies when • too ind • 
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T 

x(T) » / u dt + x(0), 



where x(T) is the position of a liquid particle at time T that passed 
through x(0) at tine 0. Figure 2.8 shows monotonlc decay of a sinuous 
free surface upstream to a nearly uniform plug flow at a distance of ' 
twice the real part of the decay rate. a RM] . downstream. 

The situation Is markedly different when capillary number Is 
decreased from 100 to I. The first two decay rates coalesce and form a 
complex pair at a wavenumber or dbout 0.5 as shown In Figure 2.3. For 
wavenumbers greater than about 7. «, 3 becomes the most slowly decaying 
rate. a p and the complex pair of decay rates, aj and a 2 , became «„ and 
tt HI» respectively. This crossover of slowest decay rates occurs only 
at finite wavenumbers. The transition from real « t to complex ^ and <x 2 
Is shown in Figures 2.10-2.12. Figure 2.10.1s a portrait of the streak- 
lines for the mode corresponding to ^ when wavenumber Is 0.1. The decay 
is purely exponential and the flow pattern Is similar to that shown In 
Figure 2.8. However, when wavenumber Is unity the slowest decaying mode, 
corresponding to oj and « 2 , Is oscillatory In the downstream coordinate, 
as shown In Figure 2.11. and the streakHnes look substantially different 
that those in Figure 2.10. At a wavenumber of 10, the mode corresponding 
to aj and a 2 Is shown In Figure 2.12 tc oscillate less noticeably in 
streamwlse coordinate than at a wavenumber of unity because the exponen- 
tial decay Is more rapid. Figures 2.13 and 2.14 show In more detail the 
flow profiles of the first two modes for selected values of wavenumber. 
Figure 2.13 shows the three components of velocity, scaled so thaj maxi- 
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mum values are unity, for two wavenumbers . 0.1 and 0.398. which are 
below the limit at which the two decay rates coalesce. Each of the 
curves Is nearly parabolic In form. 

The situation Is drastically different at higher values of wave- 
number. In Figure 2.14 the Imaginary parts of the three velocity compo- 
nent profiles are shown for wavenumbers of 0.501, 1, 2.00. 3.16, and 10, 
all above the limit at which the first two decay rates form a complex 
conjugate pair. These velocity profiles would occur at planes perpen- 
dicular to Uie x-axIs and spaced at a distance of 2it/aj apart. Tim pro- 
files become Increasingly complicated as wavenumber is Increased. When 
wavenumber Is 10, the motion Is confined to a thin layer near the top of 
the liquid film of which the thickness Is about one half of the trans- 
verse wavelength. »/10. The spatially oscillatory perturbation 1s unable 
to reach the web because viscous diffusion Is slight. Lamb's (1945) 
analysis of two-dimensional standing waves on an Invlscld moving stream 
1s analogous. His result shows that motion varies as e " M(1 "J r) when the 
wavelength Is short compared with stream depth. Thus, the velocity dis- 
turbance at a depth ir/K below the free surface Is only about 4S of that 
at the free surface. Stokes' (1851) analysis of the flow near an 
oscillating flat plate Is also analogous. He showed that If the wall 
velocity Is U Q cos(nt), the velocity of the liquid decays exponentially 
with distance from the wall: 

"(y.t) = U Q e" ky cos(nt - ky) . 




QOO*t 



a 'J.H9I3H 



o 
6 



<U 

ii "a 
o 

s e 

"a u 
c 

o 

O f 



BOOM 




000*1 



'1H9I3H 



76 




A '1H3I3H 



77 



where k s /n7ST . Thus, the liquid oscillations are confined to a boun- 
dary layer wh1v.h has a Wckness of about »/k (based on the « criterion 
used above). 

At lower capillary number the behavior of the first two decay rates 
with wavenumber Is even more complicated than In Figure 2.9/ as shown In 
Figure 2.15. The first two decay rates merge and form a complex conju- 
gate pair at wavenumber slightly greater than 0.1. At a wavenumber of 
about 0.5 the real part of this complex conjugate pair Increases rapidly 
with wavenumber and soon overtakes the next three decay rates, a, Is 
the dominant mode, Oj, for wavenumbers higher than about 0.6. Figure 
,2.16 Is a portrait of six streaklines In a How disturbed by the first 
mode when wavenumber Is unity. Although the decay rate is complex, the 
free surface elevation appears to decay monotonically with downstream 
distance because the half-wavelength, «/« I(Mg » 5.45, 1$ much longer than 
the distance necessary for the now to become nearly solid body transla- 
tion, which 1s a distance of about 2/a Rea , » 2.28. However, the flow 

would exhibit large standing waves if wavenumber were O.S because a T 

Imag 

5 0.3 Is about five times as shown in Figure 2.15 and therefore 
the half-wavelength, »/a Imag a l0 . 5 , is much shorter than the distance 
2/a Real 3 33 * In 9 enera1 ^wld be expected that a mode "appears" to 
be spatially oscillatory, i. e . to exhibit standing waves that are notice- 
able to the naked eye only when • lBn g/<» RMl > 1. 

The effect of decreasing teynolds number from 100 to 10 to I on the 
behavior of the first five spatial decay rates. * x through a 5 , as wave- 
number varies when capillary number Is 0.01 Is seen by comparing Figures 
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Figure 2.15. Behavior of the first five special decay rates, through 

ae« « wavemrsfter varies Mfce* • ioo and N- • 0 01 ' 

real part; , toa?1nary part. 8 
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2.15, 2.17, and 2.18. When Reynolds number .1$ 100, Figure 2.15 shows 
that the dominant mode Is exchanged at a wavcnuober of about 0.6 as 
described above. This outright exchange of dominant modes evolves Into 
a complicated. pattern of merging of real decay rates Into complex pairs 
and splitting 'again back to real ones when Reynolds number 1s decreased 
from 100 to 10 as shown In Figure 2.17. At even lower values of Reynolds 
number the complicated pattern of Figure 2.17. gives way to the simple 
behavior shown 1n Figure 2.18 for Reynolds number equal to 1. 

Figures 2.7, 2.9, 2 .iS, 2.17, and 2.18 show th*t the dominant mode 
can exhibit "observable" spatial oscillation - as defined by the cri- 
terion above that ^ g /^ i} > I - only when capillary number Is of 
order I or smaller and Reynolds number Is of order 100 or greater. We 
conclude .then that downstream standing waves on a coated film can be 
detected by the naked eye only when . operating and liquid parameters are 
so chosen that Reynolds number Is over 100 and capillary number Is below 
1. and furthermore the flow is disturbed by a sinuous transverse pertur- 
bation with a wavenumber of about unity (depending on capillary and 
Reynolds numbers), such as those generated by the ribbing Instability. 
However, use of optical techniques such as Hore topography (see Kheshgl 
1982) should significantly Increase the range of operating parameters 
for which standing waves could be seen. 

Comparison of Analyses. The slowest decay rates predicted by 
lubrication theory In Section 2.3 and the Oseen analysis In Section 2.4 
are compared In Figure 2.19. The Oseen prediction Is Indistinguishable 
from that of lubrication theory when Reynolds number 1s zero (as It Is 
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In lubrication theory) and when wavenumber is 0.1 or less. At higher 
wavenunber the results differ markedly, the full theory predicting much 
slower decay rates than those of lubrication theory. Fall's (1978) 
solution is shown 1n Figure 2.4 to be in good agreement with the general 
lubrication result of Equation (2.3.5) only when » ca /N 3 1s large, and 
therefore Fall's assertion that his result is valid at low capillary 
number is false.' 

Rate of Spatial Decay. The results in Tables 2.1-2.5 can readily 
be used to estimate downstream distance needed to achieve a certain 
degree of leveling- of a snail disturbance. Table 2.6 shows as an exam- 
ple the downstream distances needed for .disturbances of wavenumber 0, 
0.1, and I to decay by a factor of 10 when a coating device Is operating 
at a Reynolds number of 1 and a capillary' number of 10. As expected, 
long wavelength (small wavenumber) disturbances are the most dangerous 
and two-dimensional disturbances (wavenumber equal C) decay the most 
rapidly. Disturbances with wave numbers near unity, such as those gener- 
ated by the ribbing instability (cf. Chapter 6), decay by a factor of 10 
after traveling a few hundred film-thicknesses downstream, and therefore 
may be tolerable if there is a long enough run of the web downstream of 
the coating station. 

Construction of Robin Bound ary Condition. As mentioned in Section 
2.1. the Oseen -analysis of Section 2.4 can be used to construct a Robin 
boundary condition, i.e. a boundary condition of the third kind, for the 
outflow boundary in a three-dimensional finite element stability analysis 
of a coating How. The equations that describe the stability of a two- 
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Table 2.6. Ownstream distances for disturbances to decay by a factor 
- 10 when Reynolds number 1s I and eapflli^nSblr 5 U. 



Vavenumber, M 
(film- thickness* 1 ) 


Wavelength, fcr/N 
(film- thickness) 


Oec^y distance, ln(10)/a n 
(fUra-thickness) R 


0 


• (2-0) 


3 ' 


0.1 


62.8 


700,000 


I 


6.3 


190 



dimensional coating f!ow to three-dimensional perturbations are derived 
in Chapter 5. Equation 5.2.25, the Galerfclr. «*k form of the momentum 
equation, contains an Integral over the outflow boundary of the forming 
zone: 



L 1' I" ' (♦ J (5.n)0(Nz)] ildndz 
<0 OUT - 3n . 



(2.6.1) 



where n 1s the outward pointing normal to the outflow bound**, T' Is 
the change In the stress tensor due to a perturbation 1n the vebclty 
aid pressure field. ♦ J f«.n) 1s the j-th finite element basis function as 
defined In Section 4.2). O Is a tensor of Fourier basis functions as 
defined In Equation (5.2.10), and 6^ represents the outflow boundary 
of the fanning zone. When the flow field at the outflow boundary Is well 
approximated by the dominant mode of the Oseen analysis and n is chosen 
to be X for convenience, the stress tensor is approximated by:' 



1 • T' s A • u* . 



(2.6.2) 
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The momentum .transfer coefficient tensor Is found from Equations (2.4.9) 
aid (2.4.10) to be 

A » {1[23u(N,a 1 )/ax - p(N,o I )]/u(N,a I ) 

* JPu(N.e 1 )/»y]/u{!!. Jl ) + ^Ca u (N 1 a I )/az]MN l a I )}1 > 
+ JiC3v(M.aj)/ax]/v(M,a I ) ♦ kkfjwfN.aj )/3xJ/w(M,Bj ) 
3 {" * P 1 f.« I ^(o I )/u[ J (N, ai )s J (« I )] 
+ JCu^ (H, 0l 4 ( tt j jSjfaj J/u^tH.aj JS^taj ) J 

- k NtantHzJ^-ajm + JJ > kk) . (2>6>3) 

The smaller the real part of the dominant decay rate Oj as compared with 
the real part of the second decay rate «„, the better Is the approxima- 
tion In (2.6.2). Fortunately, the real part of the dominant decay rate 
is smaller than the real part of the second decay rate by an order of 
magnitude or more for a wide range of the parameters fn Tables 2.1-2.5. 
This Is especially true when capillary number Is greater than 0.1 as it 
Is at the onset of the ribbing instability (see Chapter 6). Equation 
(2.6.2) Is Imposed as a natural boundary condition by replacing n • T' 

mm 

1n '2.6,1) by A • u' as explained further 1n Chapter 5. Outflow boundary 
conditio.? (2.5.2) can also be used In three-dimensional finite element 
analysis of ste'acjy flows. 

2.7. Extensions of the Oseen Analysis 

The methods used 1n this chapter to construct asymptotic solutions 
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in a downstream free surface zone apply equally well to an upstream free 
surface zone or to upstream or downstream channel zones. The analysis 
of an upstream free surface zone 1s Identical to that culminating In 
Equations (2.4.34)-(2.4.52) except that the relevant roots of the char- 
acteristic equation (2.4.52) are ones with negative real part. The anal- 
ysis of the spatial decay of three-dimensional disturbances In a channel 
zone with upstream or downstream distance differs from that in a free 
surface zone only by the boundary conditions at the upper surface. The 
methods of Section 2.4 therefore can be applied to extend Wilson's (1969) 
analysis of the spatial decay of flow disturbances In a channel to in- 
clude three-dimensional disturbances. The extension is straightforward 
when inertia fs negligible (M Rft - 0). Wilson's Equation (3.3) applies, 
but with 2a Q replaced by 8 s /a* - N 2 . (Here the solid boundaries are 
at y * 0 and I Instead of y » t l as Wilson uses, which accounts for the 
factor of 2). The characteristic equation 1s then 

sln(s) - t b . • C2.7.1) 

The first root of Equation (2.7.1) 1s 

a l * H - (2.7.2) 

and higher roots can be found from Wilsons Table 1 by using the rela- 
tionship 

°i*l a /{2a oi^ + n2 • * a 1.2,... * (2.7.3) 

where a Qj is Wilson's i-th elgenvaue- a Q . Solutions at finite Reynolds 
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number could be constructed as well, but would require computer-based 
methods of analysis because of the parabolic baseflow, U(y). The Oseen 
equations for this case are 

SeW S7 + v lijJ 8 - *P + '* Z u (2.7.4) 

7 -i e0 (2.7.5) 
and the no-slip and no-penetration boundary conditions are: 

u -0 • y -0,1 . (2#7-6 , 

In this case, the methods of Section 2.4 lead to a fourth-order ordinary 
differential euqatlon with non-constant coefficients for which solutions 
are best constructed by numerical means. 
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CHAPTER 3 — LUBRICATION ANALYSIS OF THE STABILITY nc q nT vuicc 



3.1. Introduction 

The stability of steady, viscous, near.y rectilinear Hows can be 
studied by making a set of approximations that are collectively known as 
the lubrication approximation. Osborn Reynolds (1886) first proposed 
the Ideas of lubrication theory, which he used to explain some experi- 
ments on lubrication In Journal bearings done by Mr. Beauchamp Tower. 
However, it was not until 1960 that loss of stability of thin film flows 
was first Investigated in the context of lubrication theory. Pearson. 
, motivated In part by the observation of so-called brush marks forced by 
a brush spreading paint, analyzed the onset of ribbing in a wedge-shaped 
spreader. I.e. a gently sloped wedge-shaped die translating perpendicular 
to the direction of its translational symmetry and parallel to a plane 
surface. Pearson found that the results of stability analysis under the 
lubrication approximation agreed qualitatively with his experimental 
findings. Pitts and Greiller (1961) applied Pearson's ideas to estimate 
the conditions for onset of ribbing in film splitting between two sym- 
metric counter-rotating cylinders. Film splitting differs from spreading 
«n that the meniscus has no point of attachment to a solid surface. I.e. 
has no contact line, in film sp ,itt1ng whereas spreading always has a 
static concact line. As a result, the meniscus is less constrained 1n 
film splitting which causes the flow to be less stable, as compared with 
spreading. Taylor (1963) perfonned a series of experiments on cavita- 
tion in lubrication from which he deduced the Importance of capillary 
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number to the separation conditions rsar a static contact line where the 
premises of lubrication theory are Invalid. Coyne and Elrod (1971), 
apparently inspired by Taylor's work, constructed an approximate set of 
meniscus conditions based on an ad hoc model. Later, Savage (1977) 
showed how Coyne and El rod's approximate boundary conditions could be 
Incorporated into a lubrication staMU*;- analysis in his work on cylin- 
drical spreaders.' 

The chief drawback of lubrication theory when applied to partially 
submerged bearings, Ke. bearings in which separating contact lines occur, 
or to coating devices is that the boundary conditions that should be 
imposed at a separating contact line. are, not well established. Several 
boundary conditions hive been postulated, but which, if aoy of them, 
Is best remains uncertain. One of the aims of this chapter is to eval- 
uate some of the postulated boundary conditions by calculating the sta- 
bility results that they predict. The comparisons drawn here are for 
exponentially-shaped slot and knife coaters. 

The lubrication equation and postulated boundary conditions that 
apply to slot and knife coating devices are laid out . in Section 3.2. 
Integration of the lubrication equation for translatlonally symmetric 
flows is outlined in Section 3.3. Stability of the two-dimensional 
flows to three-dimensional ribbing disturbances 1s the topic of Section 
3.4. Finally, the results are discussed in Section 3.5. 

3.2. Lubrication Equation and Boundary Conditions 

Derivation of the lubrication equation for the standard case of a 
flow confined between two solid su^aces Is analogous to the derivation 
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given in Section 2.2 fcr the case where one bounding surface is a 
meniscus. Equations (2.2.11)-(2.2.16) apply here, with the exception 
that the shear stress condition (2.2.14) is replaced by a no-slip con- 
dition at the solid die: 

kl'SL 0 (3.2.1) 

Here u again represents velocity and the subscript II Indicates that only 
streamwise (x) and transverse (z) velocity components are represented. 
Following the same procedure as outlined In Section 2.2 to eliminate 
velocity In favor of pressure leads- to the standard lubrication equation 
(see e.g. Cameron 1976): 

* n -(h 3 v IlP /l2) » Oh/3x)/2 (3.2.2) 

Equation (3.2.2) differs from the result in (2.2.19) by a factor of one- 
half on the left side of the equality ihlch is due to the difference In 
boundary conditions at y * L Here h refers to the local gap thickness 
between the stationary die and the moving web and p is the pressure. 

Ideally, coating Hows are transversely symmetric. When this fs 
so, Equation (3.2.2) reduces to 

4_ (ill 4L) . I «L ' u 

dx dx\ 2 dx (3.2.3) 

Equation (3.2.3) apparently requires two boundary conditions to .specify 
a unique solution. Usually, one of these is at the -Inflow boundary. 
However, when the position of the downstream boundary, i.e. the separat- 
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Ing contact line (see Figure 3.1), is unknown an additional coition is 
needed there. In general two conditions must relate the contact line 
position, the pressure there, and the pressure gradient there. When the 
contact line is pinned at a geometrical or compositional discontinuity, 
the additional boundary condition 1s typically that the pressure at the 
known contact line location Is ambient minus capillary pressure. When 
the contact line" Is free to move along the die, the condition specifying 
contact line position is normally replaced by a condition on pressure 
gradient at the contact line. Models to determine capillary pressure 
and pressure gradient at the contact line are thus needed. 

Capillary Pressure Models. The simplest way to account for capil- 
lary pressure is to approximate the meniscus profile as an arc of. circle 
near the contact line and as a line further downstream. The contact 
angle is a parameter. The splice between circle and line is chosen so 
that free surface elevation and slope are continuous, as Illustrated in 
Figure 3.1. The dimenslonless pressure Just upstream of the contact 
line 1s then 

P . H R * (3.2.4) 

C3 

where M ca = uu/o 1s the capillary number and R Is the radius of the cyl- 
inder. The ambient pressure is absorbed Into p as i datum. Equation 
(3.2.4) 1s an approximation to the normal stress boundary condition, as 
given in Equation (1.2.5). Normal viscous stresses are always considered 
negligible in the lubrication approximation, even though they may not in 
fact be negligible In the vicinity of ths contact line where the flow 
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may be far from rectilinear. A more elaborate How model Is needed if 
the effects of normal viscous stress are to be accounted for. 

Coyne and El rod (1971) proposed an ad hoc flow model to estimate 
normal viscous stress as well as capillary pressure at a separating con- 
tact line. They assumed that contact angle Is 90° and the velocity along 
spines perpendicular to the meniscus Is a quadratic function of distance 
from the meniscus to the web. They calculated the ratios d/h. and R/d 
as functions of capillary number, where d is the gap thickness h at the 
separating contact line. Their ratio R/d can be used in place of the 
arc-of^-cyllnder approximation in (3.2.4). Salto and ScHven (1981) 
showed Coyne and El rod's results to be In surprisingly close agreement 
with their finite element results, thus lending credence to Coyne and 
El rod's ad hoc model. 

Unfortunately, Coyne and El rod's boundary conditions are not easily 
used when the separating contact angle is other than 0° or 90", the two 
cases for which they tabulate outcomes of their calculations, and espe- 
cially when contact angle is greater than 90° because then an extrapola- 
tion of their results would be necessary. An alternative to using Coyne 
and El rod's results Is to calculate the ratios d/h, and R/d by construct- 
ing finite element solutions of the exact governing equations as explained 
In Chapter 4. By so doing, Coyne and El rod-like boundary conditions can 
be generated at' all capillary numbers and contact angles. 

Pressure Gr adient Models . Several conditions that specify pressure 
gradient at the separating contact line have been proposed. Host of 
then stem from ad hoc arguments and so are at- best of dubious validity. 
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They are described briefly 1n the following paragraphs. 

Reynolds (1386) originally proposed that the pressure gradient at 
the separating contact line should be zero, I.e. 

ap/3x * 0 9 x = c . (3s2 .5) 

Here c stands for the location of the separating contact line. His 
rationale was based on the heuristic argument that the meniscus cannot 
withstand a pressure gradient for it has no force of its own by which to 
resist (this argument negfects the action of surface tension). Later, 
Swift (1932) and Stieber (1933) working Independently and employing dlf- 
. ferent reasoning arrived at the same boundary condition as had Reynolds 
almost fifty years earlier. Reynolds 1 boundary condition predicts a 
•pressure distribution such as that In. Figure 3.2a when the die is a 
cylinder, I.e. It predicts that the separation line occurs at the point 
of Incipient pressure minimum. 

Prandtl (1937) and later Hopkins (1957) proposed that the, separa- 
tion line should be at the point of Incipient reverse How, This argu- 
ment leads to the following boundary condition: 

dp/dx * 2/h 2 . Q x * c . '(3.2.6) 

This boundary condition Implies that a subambient pressure zone occurs 
Just upstream of the separating contact line, as shown In Figure 3.2b 
for the case of a cylindrical spreader. 

Taylor (1963) was perhaps the first to realize that capillary number, 
which measures the relative importance of normal viscous stress at the 




Q. 
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meniscus to capillary pressure, plays a significant role In determining 
not only the piessr-* but also the pressure gradient at the separating 
contact line, Althut-ji he did not formulate his own boundary condition, 
his paper appears to have motivated Coyne and El rod's (1969) approximate 
analysis to account for the effects of normal viscous stress* 

As already mentioned, Coyne and El rod estimated the ratios d/h w and 
R/d as functions of capillary number. The first of these ratios is a 
relationship between pressure gradient and contact line position. This 
Is easily seen by Integrating the lubrication equation (3.2.3) once to 
get: . 

- h 3 (dp/dx)/12 + h/2 * q (3.2/7) 

The first term Is the dlmenslonless flow rate per unit channel width due 
to Polseullle. or pressure-driven, now; the second term 1s the dimension- 
less How rate per unit channel width due to Couette, or shear-driven, 
flow. Therefore, an overall mass balance shows that q must be the final 
film thickness h w . Hence the pressure gradient at the contact line 1s 
given by: 



dp/dx = 6(1 - 2h /h)/h 2 



(3.2.8) 



This boundary condition predicts pressure fields similar to those pre- 
dicted by the Prandtl -Hopkins condition, as Illustrated 1n Figure 3.2b. 
The main distinction between the Coyne-Elrod and Prandtl -Hopkins condi- 
tions 1s that the former predicts a zone of recirculation just upstream 
of the separating contact line when the contact angle is 90°, the latter 
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does not. 

A cruder boundary condition thai* the three described above 1s due 
to Gumbel (1921). but often bears the name half-Sonsnerfeld. Gumbel's 
condition states that the separating contact line Is at the point of 
minimum gap width, and thus it applies solely to die configurations In 
which there is a local minimum clearance. It Is used chiefly to esti- 
mate the load capacity of journal bearings. The pressure profile that 
Gumbel's condition predicts Is illustrated in Figure 3.2c for the case 
of a cylindrical die^ 

Savage (1977)- and Dowson and Taylor (1979) have reviewed each of 
the foregoing boundary conditions. They conclude that Coyne and El rod's 
boundary condition agrees best with a' collection of .experimental data of 
the pressure distribution In bearings. However, Savage points out that 
Coyne and El rod's boundary condition is In poor agreement with Smith's 
(1975) experimental data on the pressure distribution In a lightly 
loaded Journal bearing, and especially so when capillary number is small. 

It remains to be seen which of the boundary conditions performs 
best in lubrication stability analysis. Savage (1977) chose the Coyne- 
El rod boundary condition in his stability analysis of a cylindrical 
spreader and found reasonable agreement with experimental data. However, 
no side-by-slde evaluations of the various postulated conditions have 
been made. Here, three of the boundary conditions on pressure gradient, 
namely Reynolds', Prandtl -Hopkins' . and Coyne-El rod's conditions, are 
compared for slot and knife coating when the die shape is exponential . 
In Chapter 6, the lubrication ssabillty results are evaluated by com- 
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paring with finite element stability results. 

Inflow Conditions. Three linear inflow conditions can be Imposed 
at an Inflow boundary. Of these, the first Is most commonly used. The 
three choices are described In the following paragraphs. 

The first Is to specify pressure at an Inflow boundary: 

pSp I 9 x ' b (3.2.9) 

This boundary condition is typically used In problems where the flow Is 
not premetered. 

~.e second alternative is to specify the pressure gradient at the 
Inflow boundary, which is equivalent to specifying the flow rate (see 
Equation (3.2.8)): 



dp/dx » pj 9 x » b 



(3.2.10) 



This boundary condition Is useful when the flow Is premetered, as it Is 
In many coating devices. 

The third alternative Is to .specify pressure gradient as being pro- 
portional to the difference between an external pressure and the pressure 
at the inflow boundary, i.e. the Robin condition 

dp/dx »« ( Pq . p) . . a2>11) 

Here a is a constant of proportionality and p 0 is the external pressure. 
Equation (3.2.11) is analogous to the heat flux boundary condition used 
In heat conduction problems which presumes that the heat flux at a 
solid/gas boundary Is proportional to the difference bet* en boundary 
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and ambient temperatures, the constant of proportionality being the heat 
transfer coefficient. Here the proportionality constant a Is In effec? 
a momentum transfer coefficient. 

In the following analysis pressure Is specified at the inflow bound- 
ary. A comparison given in Section 6.5 of lubrication stability analy- 
sis and finite element stability analysis is for pressure gradient spec- 
ified asymptotically far upstream. The Robin condition (3.2.11), while 
being perhaps the most physically realistic of the three alternatives, is 
not considered further here. However, the pressure and pressure gradient 
inflow conditions can be thought as' limiting cases of the Robin inflow 
condition as a Increases without bound and goes to zero, respectively. 

3.3. Two-Olmens1ona1 Base Flows 

Two-dimensional solutions of the lubrication equations are found 
by integrating equation (3.2.8). The standard solution Is (see e.g. 
Cameron 1976): 

x 

/ 
b 



P ■/ [6(1 - 2h./h)/h 2 ]dx + P0>1 (3.3.1) 



Equation. (3.3.1) contains two unknown constants which have to be deter- 
mined, i.e. the final film thickness h., and the inflow pressure p(b). 
Furthermore, the position of the contact line c is unknown. There are 
determined by choosing an inflow condition and two conditions at the 
separating contact line, as explained In Section 3.2. 

Here the case of an exponentially diverging channel is considered, 



i.e. 
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h * e x/S 

" e (3.3.2) 

The contant 1/B determines the characteristic slope In the channel and 
the inflow boundary is chosen to be at x * 0. For this die configura- 
tion. Equation (3.3.1) gives: 

p - 38(1 - e- 2x/B ) - 4Bh„(l - e- 3x/B ) ♦ p ( b ) . ( 3 . 3 .3) 

If the inflow pressure is given by (3.2.9) then p(b) Is known; if on the 
other hand inflow pressure gradient is given by (3.2.10) then h. is 
determined by (3.2.8). The contact41ne position c and either the inflow 
pressure p(b) or the final film thickness h, are specified by choosing 
two conditions at the separating contact line. 
The Reynolds condition (3.2.5) implies that 

h - * J h(c) • (3.3.4) 

as seen from the integrated lubrication equation in (3.2.8). Likewise, 
Prandtl-Hopklns' condition (3.2.6) implies that 

h -"j h(c) ' . (3.3.5) 

These conditions are Imposed here In conjunction with the arc-of-clrcle 
approximation to the meniscus profile at the contact line. For ease of 
comparison with the Coyne-El rod boundary condition, contact angle is 
chosen to be 90°. The pressure at the contact line Is given by (3.2.4) 
and the radius of curvature R Is h(c) - h^ so that 
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h(c) - h w 



Q x » c 



(3.3.6) 



Coyne and El rod's boundary condition is Imposed similarly to those 
of Reynolds and Prandtl -Hopkins. Coyne and Elrod's ratios 



tabulated as functions of capillary nuaber, 4re used to eliminate the 
final film thickness, from (3.3.3) and R from boundary condition 
(3.2.4). 

The contact line positions predicted by the three alternative bound- 
ary conditions are compared in Figure Z.Z for the case where the Inflow 
pressure is specified to be zero and the, characteristic slope 1/B 1s 0,1. 
Reynolds 4 and Prandtl-Hopklns 1 conditions predict similar separating 
contact line positions, and especially so-*hen capillary number 1s small. 
The Coyne-Elrod boundary condition predicts that the separating contact 
line should He significantly further downstream than the other two con- 
ditions, except at capillary numbers near unity where Coyne-Elrod and 
Prandtl-Hopklns conditions agree closely. 

3.4. Stability of Base Flows 

Pearson (1960) estimated the stability of a lubrication flow under 
a wedge-shaped* spreader by: (1) assuming a small, disturbance to the 
pressure Meld and to the meniscus shape; (2) using a Fourier normal mode 
representation so that disturbances of Individual wavenumber could be 
singled out; and (3) evaluating stability by means of an integral mass 



* = Vh(c) ; . b s R/h(c) 



(3.3.7) 
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balance around the separation, line, i.e. relating the rate of change In 
contact line position to the difference of Inflow rate just upstream of 
the contact line and outflow rate just downstream of the contact line. 
The method used here is nearly Identical to Pearson's. However, the 
analysis here'ls carried out by using the same three boundary conditions 
as in Section 3.3., whereas Pearson resorted to experimental data to 
determine the position of the contact line. 

If the presure field and contact line position are perturbed 
slightly they will either return tr- their undisturbed states or continue 
to move away from them. If they return the flow is said to be stable; 
otherwise it is said to be unstable. Here, the stability of a base flow 
to disturbancss of the form 

p' - c*(x) cos(Nz) (3.4.1) 
C -cco,(Ns) - {3>4>2 . 

is considered where p' and c' are snail perturbations to the pressure 
field and separating contact line position so that 

>° "'*>'. (3.4.3) 
E = C+C '- (3.4.4) 

p and 5 represent the perturbed pressure field and separating contact 
line position, respectively. A snail disturbance to the pressure field 
satisfies the homogeneous form of the lubrication equation (3.2.2): 



v 



II-< h V>-° (3.4.S) 
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Fop p' given by (3.4.1), this reduces to an ordinary differentia! equa- 
tion In 

d 2 »/dx 2 + (3h x /h)d*/dx - M 2 * - 0 (3>4>6) 

Finally, replacing h by the exponential form given In (3.3:2) results In 

d */dx 2 + (3/8) d*/dx - N 2 * - 0 . (3t4>7) 

At the Inflow boundary. * mist be zero In order to satisfy the condition 
that pressure be zero there. The boundary condition at the separating 
contact line depends on how the curvature there changes with the param- 
, eter e. 

A model Is needed to .determine how the disturbance to the contact 
line affects meniscus curvature. Figure 3.4 Illustrates the two types 
of meniscus disturbances that are considered here. The first Is a vari- 
cose deformation 1n which the radius of the cylinder Is perturbed by a 
sinusoidal function of the axial coordinate. The second is a. sinuous 
deformation In which the cylindrical axis Is perturbed as a sinusoidal 
function of the axial coordinate. The change In curvature due to the 
superposition of these two meniscus deformations when contact angle Is 
90° is 

• (3.4.5) 

Here H' is the change In mean curvature at the separating contact line 
due to a disturbance given by (3.4.2) when the relative amplitude of 
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Figure 3.4. Illustration of varicose and sinuous deformations to a 

cylindrical meniscus due to three-dimensional disturbance. 
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varicose to sinuous deformation Is given by r. The pressure disturbance 
p' must therefore satisfy the boundary condition: 



V ■ « cos(Nz)[- N 2 + i (_£_)]/H ca . c C0S(N2) Md t x a 



c 

(3.4.9) 



The last term is due to the displacement of the contact line through the 
pressure field of the base How. The boundary condition on ♦ (*) 1 $ : 

♦ * (- n 2 + ! r_L.ii/M _ Mc) a 

L R r+1 ca ~dx" 9 xac ' I3A.IQ) 

The solution of (3.4.7) that goes to zero at the upstream boundary and 
satisfies (3.4.10) 1 S 

♦ = ([- N 2 ♦ i (JL JJ/Nca - iEi£»j e Vc-x) ^V> 

R r+1 ca J, simile) ' ■"•♦•"J 

where A. = - 3/2B.and A» = /A 2 * N 2 . 

The ratio of varicose to sinuous meniscus deformation has yet to be 
specified. The lubrication approximations breaks down near the separating 
contact line because of the steep slopes there, and so lubrication theory 
cannot be used to estimate the meniscus profile there. Pearson (1960) 
was the first to realize the uncertainty In the change of meniscus cur- 
vature due to a contact line disturbance. For lack of information, he 
was forced to make an ad hoc choice of how the meniscus deforas. He com- 
pared two choices and found that his stability results were nearly the 
same for both. His first choice was to allow only the sinuous deforma- 
tion, which is to set r = 0. The second was to so choose the relative 
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amplitude of sinuous to varicose deformation that the meniscus profile 
In every x-y plane Is geometrically similar, which is to set 

r+1 B .* (3.4.12} 

Sivage's (197?) choice was also equivalent to (3.4.12). Equation (3.4.12) 
Is used In this analysis as well. 

The final step in stability analysis in the lubrication approxima- 
tion 1s to devise a criterion for Judging whether a How is stable or 
unstable to a specific disturbance. The criterion used here 1s Iden- 
tical to the one proposed by Pearson (1960). It Is that a change fn 
location of the contact line can only be due to a net accumulation (or 
loss) of mass near the contact line, which In turn must be caused by an 
Increase (or decrease) 1n inflow rate Just upstream of the meniscus 
minus the outflow rate just downstream of the meniscus. In other words, 
the criterion Is based on an Integral mass balance around the s-paratlng 
'ontact line, which is 

— fr'h(r)l - - " 3(c> 3 P* h« 

dt tc nicjj -_-£_ .„_ . (3.4.13) 

The term to the left of the equal sign represents the accumulation of 
mass on account of motion of the contact line; the terms to the right of 
the equal sign .represent the respective changes In Inflow and outflow 
rates due to the perturbation. Here 1s the change In local final film 
thickness due to the disturbance, and Is determined by the meniscus 
deformation: 



h; * e cos(fiz)[h(c)/B + r/(m)J 
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(3.4.14) 



The first term in the square brackets represents the increase- in gap 
width at the contact line because of the disturbance; the second term 
represents the decrease In radius of curvature In the x-y plane because 
of the disturbance. Small disturbances from the steady base flow either 
grow or decay exponentially 1n time, 

c = c o eU • (3.4.15) 

Combining (3.4.13)-(3.4.15) and (3.4.2) gives the following expression 
for exponential growth rate: 

x ■ - JLisLA . [h (c)/B ♦ r(m)]/h(c) (3.4.16) 

When exponential growth rate X is greater than zero the base flew Is 
unstable; otherwise ft Is stable. 

Figure 3.5 shows the stability results for exponential slot or 
knife coating as predicted by (3.4.16) when Reynolds, Prandtl -Hopkins 
and Coyne-El rod boundary conditions are used. The results show how 
growth rate varies with wavenumber at fpur values of capillary. .number. 
The results predicted by the three boundary conditions an* strikingly 
different. Each predicts that the fastest growing disturbance occurs at 
a nonzero wavenumber, but the maximum value of growth rate predicted at 
a given capillary number varies widely with meniscus condition.' The 
most important distinction is that the Coyne-Elrod condition predicts 
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that the flow 1s unstable only at the highest capillary number while the 
others predict Instability eyen at the smallest capillary number. 

3.5. Discussion 

Several researchers, most notably Pearson (1960). Pitts and Greiller 
(1961), and Savage (1977), have shown that lubrication stability predic- 
tions are In some cases In qualitative agreement with experimental 
results and they have concluded that the mechanism of ribbing can be 
explained by lubrication theory. However, there are several uncertain- 
ties In the lubrication model, not the least of which are the appropriate 
boundary conditions to choose at the separating contact line, which 
strongly affect the stability predictions. As shown In Chapter 6, the 
uncertainties are most likely due to failure of the lubrication assump- 
tions in the vicinity of contact line where the flow Is generally far 
from rectilinear. These uncertainties force the practitioner to resort 
to guesses as to the conditions that apply there. The uncertainties 
grow particularly severe when a three-dimensional disturbance to a base 
flow is considered because little or nothing can be determined about the 
manner in which the meniscus Is deformed by lubrication theory alone. 

The apparent deficiencies of stability analysis 1n the lubrication 
approximation can perhaps best be overcome by solving the exact equations 
of motion, as given In Section 1.2, and then analyzing the stability of 
those solutions. This Is done by the Galerkln/flnlte element method in 
Chapters 4 and 6. The results of the finite element stability analysis 
of exponential slot and knife coating reported In Chapter 6 are strikingly 
different than those of lubrication theory using Coyne and El rod '.s bound- 
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ary condition, as shown In Figure 6.12. One oust conclude therefore that 
lubrication stability results are unreliable, and especially so when, 
expeplnantal or other data 1s unavailable to guide In the selection of 
conditions at a separating contact line. 
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CHAPTER 4 - A HALYSIS OF STEADY TWO-OIMENSIONAL SLOT AND KNIFE COATING 
FLOWS ! : — 

4.1. Introduction 

Slot and knife coating devices are among the most commonly used by 
coatings Industries to apply thin liquid layers onto moving webs, i.e. 
flexible sheets. A typical slot or knife coater consists of a solid die 
which may be parallel or inclined to a moving web supported on the oppo- 
site side by a back-up roll, or 1n some cases stretched between two rolls 
not far apart. Figure 4.1 shows a knife coating device In which the die 
and web form a channel which has both parallel and diverging sections. 
This chapter is devoted to finding steady two-dimensional solutions to 
the equations that govern slot and knife coating Hows, as laid out In 
Chapter I. The stability df these flows to two- and three-dimensional 
disturbances is investigated 1n Chapter 6. 

The flows considered here are of Newtonian liquids with uniform 
properties. The physical properties of the liquid and the operating 
parameters can be conveniently expressed in terms of five dlmerislonless 
parameters: (1) Reynolds number. Re s Ud/v, Indicates the relative 
Importance of inertia as compared with viscous forces; (2) flow rate. 

Q = h/d, is essentially a draw-down ratio; (3) Stokes number, St = 
2 

gd /vu, expresses the relative magnitude of gravity as compared with 
viscous forces; (4) capillary number, Ca s uU/a, 1s the ratio of a 
characteristic normal viscous stress to a characteristic capillary pres- 
sure: and (S) the three-phase contact angle that the liquid makes with 
the upper solid die, or equivalent^, the position of the contact line. 
The first three of these parameters are associated with the flow field; 
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the last two are associated with the Interface- Here, U 1s web speed, 
d Is upstream gap thickness; v Is kinematic viscosity; Is the liquid 
film thickness far downstream from the slot or knife coater; g is the 
acceleration of gravity; u Is <tynam1c viscosity; and a Is .surface ten- 
sion. The effects of gravity are usually inconsequential In slot and 
knife coating and sc Stokes number 1s usually much smaller than unity. 
Here It Is taken to be negligibly small. 

Figure 4.1 shows that slot and knife coating flows, like many other 
coating Hows, contain thrse -are or less distinct flow zones. In the 
upstream channel zone the flow departs only asymptotically from recti- 
linear How, so that th- governing equations are the momentum and mass 
conservation equations' linearized about the rectilinear regime that pre- 
vails far upstream. Likewise, in the downstream free surface zone the 
flow departs only slightly from rectilinear motion -1n this case plug 
flow, I.e. solid body translation - and so the linearized form of the 
governing equations apply there as well. However, In the Intermediate 
region, here called the forming zone, the flow field rapidly evolves 
from channel to free surface How and the full nonlinear form of the 
governing equations must be retained. 

The strategy for analyzing the steady two-dimensional How field in 
slot or knife coating is illustrated in Figure 4.2. Asymptotic solutions 
that describe the How fields 1n the channel and free surface zones were 
constructed by Wilson (1S59) arsd Hlggins {1982). respectively both by 
means of perturbation analysis. How-fields 1n the forming zone are 
constructed by Galeriin/finlte element analysis,. a method which relies 
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on the modern digital computer to overcome what would otherwise be an 
Insurmountable algebraic problem. Earlier theses by Orr (1976) and 
SilHman (1979) demonstrate that the Galerkin/finite element method is 
ideally suited to handle. viscous Free surface flow problems with irregu- 
lar boundaries; this is further documented in articles by Sllliman and 
Scriven (1980), Saito and Scriven (1981), and Kistler and Scriven (1982). 

Two important considerations are the following: (1) where should 
the inflow and outflow boundaries of the forming zone be placed? (2) 
What boundary conditions should be Imposed there? The first question 
arises because open-flow boundary .conditions are known only asymptoti- 
cally, yet when a computer-based method such as Gal eric in/ finite elements 
is used the inflow and outflow boundaries oust be located at finite 
rather than asymptotically large distances. The second question is one 
of choosing the boundary conditions that best approximate the interac- 
tions that the forming zone has with the surrounding zonesJ The choices 
considered here are the three general types of linear boundary conditions: 
(1) a Oirichlet condition, i.e. setting the dependent variables, would 
require that the asymptotic velocity profile already be attained at a 
finite location; (2) a Neumann condition, i.e. setting gradients of the 
dependent variables, would require that some final How state (not nec- 
essarily the actual one) be attained at a finite location; and (3) a 
Robin boundary condition, I.e. choosing a proportionality between depen- 
dent variables and their gradients, would require only that at a, finite 
location the now decay towards the final regime at the asymptotic rate. 
The issues of choice and location of open-How boundary conditions are 
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investigated here In two examples: developing plane Couette flow; and 
film formation on a moving web downstream of a slot or knife coater. The 
results of the first example have Implications for the inflow boundary of 
slot and knife coating flows • 

In section. 4.2 the Galerkln/finite element forms of the equations 
governing viscous free surface Hows are laid out and alternative open- 
flow boundary conditions are displayed. Section 4.3 focusses on how the 
choice and location of the outflow boundary condition affects the accu- 
racy of the results in the examples of developing plane Couette flow and 
film formation on a moving web. The main findings are summarized and 
discussed 1n Section 4.4. 

4.2. Galerkin/Finite Element Formulation 

Equations (1.2.1)-(1.2.7) govern viscous free surface Hows when 
liquid properties are uniform. These equations are solved here by the 
Galerkin/finite element method. The procedure is: .(1) the unknown 
velocity, pressure, and free surface elevation are expanded In terras of 
finite element basis functions, each of which Is nonzero only over a 
small subdomaln; (2) the weak forms of the governing equations are con- 
structed by choosing the basis functions as weighting functions and 
integrating over the domain; (3) the divergence theorem is applied to 
the governing equations wherever It is advantageous to reduce the order 
of the highest derivatives, in this case from second-order to first- 
order; (4) natural boundary conditions, I.e. ones that involve deriva- 
tives of the unknown variables, are Introduced into the boundary Integral 
that results from application of the divergence theorem; and (5) essen- 
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tial boundary conditions are Imposed by replacing the equations weighted 
by trial functions associated with boundary nodes by the respective 
boundary conditions that apply there. Each of these steps Is standard 
In finite element theory (see Strang and Fix 1973); the following para- 
graphs explain them In greater detail. 

Basis Function Expansions. The unknown velocity, pressure, and free 
surface elevation are each expanded In a set of finite element basis 
functions. Each basis function Is nonzero only on a small subdcnain of 
the entire flow domain, and Is expressed as a simple polynomial of the 
Isoparametric coordinates 5 and n which act as local stand-ins for the 
global rectangular coordinates shown In Figure 4.1. Basis functions are 
defined not only by choosing the order of the polynomial by which they 
are expressed, but also by choosing the manner in which the Joaain Is 
subdivided into elements, for this latter choice determines the subdomaln 
over which each trial function Is nonzero. 

The flow field unknowns are expanded in basis functions: 

i'Hl ♦ ,{5 ' n) (4.2.1) 

P " p 1 (4.2.2) 
h - h, p'u.n) (4<2<3) 

Here u is velocity, p Is pressure, h is free surface elevation. , ^, 
and P f are finite element basis functions as specified below; u v p i§ 
and h { are unknown coefficients which determine the velocity and pres- 
sure fields and the free surface elevation in the finite element approx- 
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filiation; and (C,n) are Isoparametric coordinates which are defined on 
each element by an Isoparametric transformation: 

(x.y) * (x 0> y 0 ) t * f U.n) , (4.2.4) 

where (* 0 .y 0 ) is the location of the i'-th node. The purpose of the iso- 
parametric transformation Is to map deformed elements Into standard 
squares dS shown in Figure 4.3, which facilitates evaluation of the 
integrals below., TVie standard procedure of mixed Interpolation (see, 
e.g., Huyakorn et al.) 1s adhered to in this analysis. I.e. velocity 
basis functions are chosen to be higher order than pressure basis func- 
tions. Here, velocity basis functions are standard biquadratics and 
pressure basis functions are standard bill nears (see. Strang and Fix 
1973). Tessellation of the domain Into elements Is illustrated In 
Figure 4.2. The location of velocity and pressure nodes in an element 
Is illustrated In Figure 4.3. Free surface elevation basis functions 
are chosen to be standard Hermlte cublcs (Strang and Fix 1973). 

Galerlcln/Finlte Element Residual Equations . Residual equations of 
the momentum conservation equation are formed by weighting Equation 
(1.2.1) by the velocity trial functions, Integrating over the domain, 
and applying the divergence theorem to the stress term. The result is: 

Rj*/ (Re.V .(uuJ^U.n) ♦ T . 7* f (e, n )f d5dn 
Q - 3U,n) 

- / n-T^(C,n)|idY (4.2.S) 
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M 

Here a f Is the, 1-th momentum residual equation; T Is the stress tensor, 
and under the conditions here that the liquid Is Newtonian and has uni- 
form properties It Is given by T = Vu ♦ (Vu) T - p(_H ♦ Jj); 3(x.y}/3( 5 , n ) 
1s the Jacobian of the Isoparametric transformation; s represents arc- 
length along the Interface; Y Is a parameter representing either c or n; 
and o and sa are the flow domain and Its boundary, respectively. Resid- 
ual equations of the mass conservation equation are formed by weighting 
Equation (1.2.2) .by the pressure basis functions and integrating over 
the domain. The result is 

»r»/{(v.«)^,n)|ili4dcd„ r (4 . 2 . 6) 



0 3(C,n) 



Here ftf is the 1-th continuity residual. Area integrals in Equations 
(4.2.S) and (4.2.6) are evaluated by standard three-by-three point 
Gaussian quadrature In the standard square reference element shown in 
Figure 4.3; the boundary Integral in (4.2.5) Is evaluated by standard 
three-point Gaussian quadrature along one of the edges of the square 
reference element. 

Adherence Bo undary Conditions . Adherence boundary conditions are 
imposed by setting the velocity coefficients corresponding to nodes 
which He on solid boundaries so that local, liquid velocity is the same 
as local boundary velocity, as Illustrated In Figure 4.4. The web veloc- 
ity is - t, so - 



—i 



(4.2.7) 
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for each node on the web. likewise, the upper solid die Is stationary; 
thus 

4"* (4.2.8) 

for each node on the die. Equations (4.2.7) and (4.2.8) replace the 
momentum residuals (4.2.S) at nodes along solid boundaries. 

Free Surface Boundary Conditions. Balance of nonnal and shear 
stresses and lack of nonnal velocity component at the liquid/gas Inter- 
face Impose three boundary conditions there, as described 1n Section 1.2. 
Here, the gas phase Is approximated as being invlscld and Is assigned a 
datum pressure of zero. The shear stress balanca is imposed as a natural 
boundary condition by replacing the boundary Integral term In the momen- 
tum residual equation (4.2.S) by 

/^nnn :T ^( 5 ,n)|idC . - (4 . 2 . 9) 

Here, «Q FS Indicates the free surface portion of the boundary and n is 
the outward-pointing unit normal to the Interface. The normal stress 
balance Is imposed as an Independent equation: 



1 S ' ! inn : T * f (£.nl ♦ fn»" l f, /{ith^Wi^Ji 

{♦<(5.n )Ca" Vf^> 1/2 }1°1 (4.2.10) 



*f ' ! (nn : T ♦'(£,„) ♦ [ca" l h -/(l+h^ja^U.nJ/ax} 3x 
™FS . 3€ 



OUT 



65a OUT and <an a refar to *■* intersection of the free surface with the 
outflow boundary and to the contact line, respectively. The last term 
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of Equation (4.2.10) Is dealt with below. The third boundary condition, 
which expresses the Impenetrability of the Interface, i.e. the kinematic 
boundary condition, 1n Imposed as an essential condition: 

' (4.2.11, 

for each node on the free surface. Equation (4.2.11) replaces the y-' 
component of the momentum residual (4.2.S) for each 1 corresponding to 
a free surface node. 

Contact Line Boundary Condition. Two alternative contact line 
boundary conditions are considered here. The simpler Is to prescribe 
the contact line position as an essential condition, which corresponds 
to the situation when a meniscus is pinned at a geometric or compost- 
tlonal discontinuity. In this case the boundary condition is 

V W . (4.2.12) 

for node 1 at the contact line, h^x,) being the elevation of the upper 
solid surface there. The second alternative Is to prescribe the angle 
that the Interface makes with the upper solid die, denoted here as 8 and 
defined by the angle between the outward pointing normals to the liquid 
and the solid die at the contact line. This is done using 

^-tan[l80- + tan- l (^)-e] (4 . 2 . 13) 

to replace dh/dx in the last term of Equation (4.2.10) at the contact 
line. 



126 

Which of the two contact line boundary conditions, Equation (4.2.12) 
or (4.2.13), Is selected makes little difference in a steacjy two- 
dimensional flow calculation, for given a contact line position there 
corresponds a contact angle and vice versa. On the other hand, it 1s 
convenient to be able to specify directly a contact angle when it is the 
known quantity, rather than iteratively choosing a contact line position 
until the one is found that satisfies the desired contact angle- However, 
the most important difference between prescribing contact angle ;r.d con- 
tact line is in the stability of the resulting How, as described 1n 
Chapter. 6. 

Open-Flow Bo undary Conditions, The slot or knife coating flows that 
are considered here are mathematical Idealizations because the domain is 
taken to be unbounded in both the upstream and downstream directions. 
In reality liquid is supplied at some finite distance upstream and the 
coated liquid eventually solidifies at some finite distance downstream 
because it is chilled, polymerized, or dried, nevertheless, these dis- 
tances are often large enough that the mathematical idealizations are 
justified. However, the Galerfc1n/f1n1te element scheme relies heavily 
on the digital computer, which is limited in speed and capacity, and 
hence the unbounded domain must be shrunk back to a finite size. One 
way to view the shrunken domain, I.e. shrunken from the unbounded domain 
of the mathematical Idealization, Is illustrated in Figure 4.2. Finite 
element basis functions are used to expand the How field unknowns In 
the forming zone. Upstream and downstream of the forming zone are 
respectively the channel zone and the free surface zone. How 1n the 
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channel zone can be made as nearly rectilinear as desired by choosing 
the Inflow boundary. I.e. the boundary between channel and forming zones, 
far enough upstream. Likewise. How in the free surface zone can be made 
as nearly plug-like as desired by choosing the outflow boundary, I.e. the 
boundary between free surface and fonnlng zones, far enough downstream. 
How far Is far enough has to be tested. Furthermore, the conditions to 
be Imposed at these stand-in boundaries should represent as accurately 
as conveniently possible the Interactions of the How Interior and 
exterior to the boundaries.' The more accurately they represent the 
Internal/external now Interactions; the smaller need be the finite ele- 
ment domain. This Is the criterion that Is used to Judge among possible 
boundary conditions: that, boundary condition is "best" which allows the 
forming zone to be the smallest given. the accuracy desired. 

Each of the three types of linear boundary conditions are considered. 
The first is a Oirtchlet condition In which the asymptotic velocity pro- 
file and free surface elevation are specified at the open-flow, boundary 
nodes: 

±'!L'l (4.2.H) 
h * h - = 0 -(4.2.15) 

Here ^ and h„ represent the asymptotic velocity profile and free surface 
elevation, respectively. Equations (4.2.14) and (4.2.15) are Imposed as 
essential boundary conditions. The second is a Neumann condition In 
which the asymptotic normal traction and free surface elevation gradient 
is specified at the open-flow boundary- nodes: 
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I- [T - TJ ■ 0 { ^ 2A6) 
dh/dx = 0 (4-2>l7) 

Equations (4,2.16) and (4.2,17) are Imposed as natural boundary condi- 
tions by using them to replace the normal traction term at the outflow 
boundary in the momentum residuals (4.2.5) and free surface gradient 
term at the Intersection of the free surface and the outflow boundary Ui 
the normal stress residuals (4.2,10), respectively. An alternative to 
Equation (4.2.16) 1s to Impose the asymptotic normal velocity gradient 
at the open-flow boundary: 

i- [vu - *J ■ 0 . (4 . 2#18) 

Equation (4.2.18) would be the appropriate boundary condition If only 
the laplaclan of velocity had been Integrated by parts in forming the 
momentum residuals. The third 1s a Robin condition 1n which normal 
traction Is taken to be proportional to velocity and free surface eleva- 
tion gradient 1s taken to be proportional to free surface elevation: 

-1 ; [T - Tj = A - (a -jj (4.2.19) 

- dh/dx = cfh hj (4.2,20) 

Here A is a ttyadlc, a second-rank tensor which depends on no more than 
the fllmwlse coordinate y f and c Is a constant. A and c should be chosen 
so that velocities and tractions are as nearly continuous as possible at 
the open-flow boundaries, as indicated in Figure 4.4. Equations (4.2.19) 
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and (4.2.20) are Imposed as natural boundary conditions In the same way 
as the Neumann condition In Equations (4.2.16) and (4.2.17), An alter- 
native to Equation (4.2.19) is 

ta-VitJ-B. («,-£.] ( 4.2. 21 , 

The asymptotic analyses of Hlgglns (1982) and Wilson (1969) show ' 
how small perturbations from final (initial) How states decay to zero 
with distance downstream (upstream) for the cases of film formation on 
a moving web and developing plane rectilinear flow, respectively. The 
results of these analyses can be used to choose the dyadic A (or 8) and 
the constant C. Both asymptotic solutions have a series fora: each, 
successive term In the series represents an increasingly complicated How 
mode which has Its own decay rate. The most important mooes are the ones 
which decay most slowly with distance, because they persist furthest 
from the forming zone. When a single mode decays much more slowly than 
any of the others - and this is usually the case in film formation on a 
moving web - the dominant mode can be used to construct the dyadic A (or 
8) and the constant C, for example: " 

As 1iC23u}/3 X -pJ)/uf] ♦ JlCOuJ/3y)/uJ] ♦ Jj[(3vJ/3y)/vJ] (4.2.22) 
C , [(h*/dx)/hj] , (4<2>23) 

or if the alternative formulation is used so that only velocity gradients 
appear 1n the boundary integral: 
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B s mOu*/8x)/u*] ♦ ii[(av*/3x)/vj] 

" ( ll*il)W (4.2.24, 

Here, uj and v* are the velocity components of the dominant mode from 
the asymptotic analysis; p* -f, the pressure field of the dominant mode; 
h. ♦ h\ Is the free surface elevation of the dominant 'rode; and « Is the 
exponential spatial decay rate of the dominant mode. The quantities in 
Equations (4.2.22)-(4,2.23) have been so chosen that boundary conditions 
(4.2.19M4.2.21) would be satisfied Identically 1f the dominant mode of 
the asymptotic cnslysls exactly dascrlbed the flow field at the open-flow 
boundary. 

Wilson's (1969) analysis of developing rectilinear channel flows at 
low Reynolds number shows that there are two dominant modes that decay 
spatially at the same exponential rate, butr'they are modulated by sinu- 
soidal spatial osclMac.ons that are 90° out of phase. Therefore, both 
modes are needed to construct the dyadic A and constant C. The method 
used here is to select coefficients » x and a 2 so that the vector equation 

is satisfied locally. Then, Equation (4.2.19) is replaced by 

where ± and ^ are defined by (4.2.21) for u* and u|. respectively. The 
same method can be used when boundary condition (4.2.21) is applicable. 
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Newton's Iteration. Newton's Iteration procedure was used to solve 
the finite element residua! equations (4.2.5)-(4.2.20). Newton's method 
and Its theory are described In a number of books on numerical methods 
for solving systems of equations, among which those by Traub (1964), 
tsaacscn and Keller (1966), Oden (1972), and Ortega and Rhelnboldt (1970) 
are especially useful . The basis of this method is simply a Taylor 
expansion, of the set of residual equations, the Nth of which Is denoted 
here as R 1 as defined by: 

N'lf-'W # .(4.2.27) 

•Here Itf* and are respectively the x- and y-components of the momen- 
turn residual £»,; n Is the. number of velocity nodes: a is the number of 
pressure nodes; s Is the number of free surface elevation unknowns. If 
6' 1s the 1-th nodal unknown, as defined by 

f'lf S K W-'-'VPl V h l h s f • (4.2.28) 

then Newton's Iteration Is defined by 

j(n) ( (n*l) (nh . R (n) ' 
J1 1 1 ' J (4.2.29) 

Here, the superscript Indicates the Iteration noAar and Jjj !s the ji-ih 
entry In the Jacobian matrix, as defined by: 

,(n) jpfn)... . 
J jt - 3R j • (4.2.30) 
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The system of linear equations expressed by (4.2.29) was solved by 
using Hood's (1976) subroutine FRONT. SHIIman (1979) used the same 
Iteration procedure as (4.2.29) In his Investigation of extrusion coating. 
However, he did not. correctly evaluate the derivatives of the residual 
equations with. respect to the free surface elevation unknowns. Saito 
andScrlven (1981) later showed how those derivatives could be evaluated 
In the reference domain of the Isoparametric mapping. Saito and Scrlven 
also demonstrated that their convergence rates were asymptotically quad- 
ratic, as theory predicts, whereas Silllman's convergence rates were 
only linear due to the errors In his Jacoblan. Asynptotlc quadratic 
convergence was also attained in this work, verifying that the Jacoblan 
matrix was evaluated correctly. 

Initial approximations, e} 0 '. to start up Newton iteration were gen- 
erated in two ways: (1) velocities were set to zero everywhere except 
at the moving web where the adherence boundary condition was imposed, 
and the free surface profile was chosen to be a straight line that passed 
through the contact line and a point lying on the outflow boundary at 
y ■ h„; and (2) a previous solution was used as the Initial guess. The 
first method was used chiefly for start-up. It usually required from 5 
to 7 Newton Iterations to converge when the criterion was 

1°" 4 • (4.2.31) 

The second method was used for continuation once the first solution had 
been attained and normally required 3 to 4 Newton Iterations. Each 
Newton Iteration required only about 5.4 central processor seconds on 
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the University of Minnesota Cyber 74, and thus the computations were 
relatively Inexpensive (typically $4.00 for one solution). 

4.3. Comparison of Outflow Boundary Conditions 

Figure 4.S Illustrates the basis used here to compare the contend- 
ing boundary conditions* Four domain lengths were chosen fn'each 
situation: for the channel, lengths of 8, 6, 4, and 2 channel widths 
were chosen; for the slot or knife coater, lengths measured downstream 
of the contact line were 12, 8, 6, and 4 channel widths. In both cases 
shorter domains were formed by chopping off elements from the longer ones. 
Therefore, the shorter the domain the fewer the degrees of freedom and 
•consequently the cheaper the computations. It was found that computa- 
tion costs were nearly a linear function of the number of elements used. 
Olrlchlet Inflow conditions were used In both cases. 

The standard of comparison was always chosen to be the solution with 
the Robin condition and the longest domain. This might seem to bias the 
results In favor of the Robin condition, but 1n fact did not because the 
three types of boundary conditions generally agreed closely on the long- 
est domain. In order to' test the Influence that the boundary condition 
has on the How further upstream, comparisons vere limited to the shaded 
portion of the domains In Ffgure-4.5, four comparisons 1n all. These are 
as follows: (1) maximum difference In the streamwise velocity component, 
u; (2) maximum difference In the cross-stream velocity component, v; (3) 
maximum difference In pressure; and (4) difference in free surface ele- 
vation at the downstream edge of the shaded region (knife coater geometry 
only). 
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Figure 4.6 shows the results for developing plane Couette ^.ow when 
inertia is negligible, i.e. Reynolds number is negligibly small. All 
three boundary conditions give comparable accuracy in velocity components, 
but the DiHchlet condition gives I to" 1.5 fewer decimals of accuracy in 
pressure. Thus for developing plane Couette flow, in which disturbances 
decay rapidly, the Robin and Neumann conditions a* superior, and com- ' 
petitive with each ether. 

The ordlnates in Figure 4.6 (also Figures 4.7 and 4.8} represent 
the number of digits to the right of the decimal point which are in 
agreement with the standard, i.e. the solution with the Robin condition 
on the longest domain. The ordinate values were calculated by taking 
the negative of the base-ten logarithm- of the difference or maximum dif- 
ference of a property at corresponding node points of the two solutions 
being compared. 

In film formation on a moving web the Rdbin condition gives the best 
accuracy when capillary forces dominate, i.e. when capillary number is 
small as shown in Figure 4.7. The Increase in accuracy with the Rabin 
condition over the others .1s about 2 decimals for velocity components, 
I to 2 decimals for pressure, and 0.5 to 2 decimals for free surface 
elevation. These are significant Improvements in accuracy, and so the 
Robin condition is to be preferred in these circumstances. 

Figure 4.8 shows that when inertial forces dominate over viscous 
forces, I.e. when Reynolds number 1s appreciably greater than unity, the 
situation changes dramatically. Here the Oirichlet condition has been 
discarded from the comparison because It led to solutions with large un- 
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physical standing waves and furthermore the convergence rate of Newton's 
Iteration was slow when It was used. The other boundary conditions, 
Robin and Neumann, are campetitlve. the Robin condition having a slight 
edge. The reason that the two are close Is clear: the spatial decay, 
rate of the most slowly decaying rode decreases nearly to zero as the 
inertia effect grows large. This is because viscosity provides the sole 
mechanism for transporting momentum from the liquid to the moving web, 
and the effect of viscosity becomes relatively snail -«hen Reynolds 
mmber becomes large. As a. result, the constants of proportionality in 
the Robin condition become very smaM at large Reynolds number, causing 
It to degenerate almost all the way to a Neumann condition, which Is ,the 
limiting case. Thus, it Is not unexpected that the two results should 
be similar. 

The chief advantage of the Robin condition Is not that It improves 
accuracy at a given domain length, but that shorter domain lengths can 
be used without loss of accuracy. Tor example, (f two-decimal accuracy 
1s sought in the velocity components and free surface elevation for the 
case of film formation on a moving web. the minium domain length, needed 
when each of the three boundary conditions Is used can be compared. 
Figure 4.9 shows that, as expected, the Robin condition allows shorter 
domain lengths to be used - about a factor of 2 shorter when' inertia is 
a small to moderate effect, but decreasing to about a factor of 1 1/3 
when Inertia Is large. 

Figure 4.10 shows the performance of the three boundary conditions 
at vanishing Reynolds number as the ratio of capillary to viscous forces 



140 





Z 



as 
o 



o 



© 
o 



*J c 
U £ 

'Z .o 

> O) 

«. * 

u > 

43 O 

3 ^ 

U (0 

u 

<o c 
o 

2 c 
e o 

u -u 

01 <0 

"? i 



E C 

o5 

> 

* r- 

tn 01 
in 

4) 01 
u o 
01 <o • 
c «*-^» 
WO 

in 3 • 
£ trt O 
+* 

o oi cn 

C Ot*F- 
OJ 

c -o -O 
— C 6 
« « 3 

E c 
O vi 

^ c C 

e oi <o 

a c *— 
E 0«- 
Q.t- 
C E CL 

o <o 
3: <j u 



01 
3 



H19M31 WnWIMIW 



142 

is varied. Again, two-decimal accuracy in the velocity ecspCRsnts and 
free surface elevation is sought. In tills case the Robin condition 
allows more than a four- fold decrease in domain length compared with the 
other two conditions. This factor, however, falls rapidly when viscous 
forces rise and It levels off at about 1,25. The results 1n Figures 4.9 
and 4.10 Indicate that the Robin condition always allows shorter domain 
lengths to be used for the given accuracy specified, but that the advan- 
tage 1s greatest when capillary forces dominate over bath Inertlal and 
viscous forces. 

4.4. Discussion 

The results of Section 4.3 show that a Robin outflow condition is 
superior to Dlrlchlet or Neumann conditions in developing channel How 
and film formation on a moving web. but the degree to which it is supe- 
rior varies widely with the circumstances^ In developing channel How 
the advantage of a Robin condition over a Neumann condition 1s marginal, 
both being superior to a Dlrlchlet outflow condition. In film formation 
on a moving web the advantage of a Robin condition Is greatest when 
capillary forces dominate over viscous and inertlal forces. One major 
advantage or a Robin outflow condition is that the degree to which the 
domain needs to be shrunk or expanded to obtain maximum computational 
efficiency without loss of accuracy as parameters are changed Is modera- 
ted In comparison with the situation when a Neumann or Dlrlchlet con- 
dition is used. For example, Figure 4.10 shows that the domain length 
measured downstream of the contact line can be maintained at about 2 gap 
thicknesses when a Robin condition is used and capillary number Is varied 
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over two orders of magnitude. By comparison, the domain length would 
have to vary by about a factor of about three ff Neumann or Dirichlet 
conditions were used. 

The question of which boundary condition Is "best" at the Inflow 
boundary of slot or knife coating has not been addressed directly, but 
some conclusions can be drawn from the foregoing results. The Inflow 
boundary of slot and knife coating flows 1s analogous to the outflow 
boundary In developing channel flow, and especially so at low Reynolds 
number when momentum convection Is not Important. This suggests that 
the type of inflow boundary condition plays only a minor role 'in affect- 
ing accuracy of computations. Moreover, when Reynolds number is not 
negligible the upstream Influence of a downstream departure from recti- 
linear How Is less than it would be at negligible Reynolds number. A 
Oirlchlet Inflow boundary condition is likely "best', and especially so 
at moderate to high Reynolds numbers. 

Figure 4.11 shows the results of two How field calculations when 
contact angle Is specified. The. arrows show the relative magnitude and 
direction of local velocity. The Inflow condition is Dirlchlet and the 
outflow condition Is Robin, m neither case has the How profile reached 
solid body translation at the outflow boundary, yet the arrows represent- 
ing the local velocities near the outflow boundary appear to be smoothly 
approaching the asymptotic state. This demonstrates the power of the 
Robin outflow condition as given in Equations (4.2.19) and (4.2,20). 

Figure 4.11 also shows the sensitivity of contact line position to 
capillary number. As capillary number . 1 ncreases (here by a factor of 
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2), the contact line Invades further upstream. Furthermore, the meniscus 
curvature Increases with capillary number. These features of the flow 
are discussed further In Chapter 6 In connection with the stability of 
slot and knife coating. 

Figure 4.12 shows some streamline contours of a similar flow but on 
a longer domain. The streamlines are parallel and quadratlcally spaced 
at the upstream boundary. Indicating plane Couette flow. At the outflow 
boundary the streamlines are parallel and evenly spaced, Indicating solid 
body translation. The streamlines fan out slightly In the diverging 
portion of the channel, rapidly converge as the meniscus forms, then 
converge mere and more gradually, finally approaching parallel and evenly 
spaced asymptotes. Although this flow appears inoccuous, it is shown in 
Chapter 6 to be unstable* to a ribbing disturbance. 
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CHAPTER 5 - FINITE ELEMENT ANALYSIS OF COATIKG INSTABILITIES 
5.1. Introduction 

Previous analyses of the fluid mechanics of coating flows have been 
Incomplete. The existence of steady two-dimensional Hows has been 
demonstrated, but not their stability or uniqueness. The stability and 
uniqueness of static equilibria, gyrostatlc equilibria, and simple rec- 
tilinear flows have been addressed (see Brown 1979, Brown and Scriven 
1980a,b, and Orazin and Reid 1981, for examples), but not of complicated 
two-dimensional Hows like those typical In the deposition of a coating 
onto a moving web. These Hows have become amenable to analysis only In 
recent years by vlrture of the advent of the modern digital computer and 
sophisticated computer-based schemes for handling elliptic partial differ- 
ential equations - primarily the finite element method. How to harness 
the finite element method to undertake the stability analysis of these 
flows Is only now coming Into focus. The origins of finite element sta- 
bility analysis in fluid dynamics date back to Orr's (1976) thesis 
research; analysis began with Brown's Investigations of the shape and 
stability of static menisci as recorded in his thesis (1979) and paper 
by Brown and Scriven (1980a). The analysis then turned to shape, stabil- 
ity of rotating drops (Brown and Scriven 1980b). Brown's s]tabHity anal- 
ysis of simple static and gyrostatlc "flows" depended on an energy mini- 
mization principle, and thereby his analysis resulted in a symmetric 
eigenproblera. On the other hand, "dynamic" free surface flows usually 
have no minimization principle because they are not conservative, and 
thereby the study of their stability leads to an asynmetrlc elgenproblem. 
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Although the Ideas and methods of finite element stability analysis 
carry over from "static" or "gyrostatlc" flows to "dynamic" flows, the 
methods to perform the elgenanalysls do not. Because techniques to solve 
asymmetric elgenproblems are not well developed at present, finding means 
of elgenanalysls necessary to determine the stability of a free surfaca 
flow was one of tha major challenges of this research. 

Pioneers Infinite element stability analysis and undamped vibra- 
tions In mechanical structures were Gallagher and Padlog (1963), Martin 
(1965), Kapur and Harry (1966), Anderson et al. (1968), Carson and Newton 
(1969). Gupta (1969. 1970, 1972. 1973). Jennings and Orr (1971), Oong 
et al. (1972), and Bathe and Wilson (1972, 1973). These early analyses 
concentrated on solving a symmetric generalized elgenproblem, i.e. an 
elgenproblem of the form 

J 1j X J " XH iJ X j ' - CS.l.ll 

where and are the respective components of symmetric stiffness 
and mass matrices, Xj is the J-th component of an eigenvector, and \ Is 
the corresponding eigenvalue. The stiffness matrix, also called the 
Jacobian matrix, Is an array of sensitivities of residuals to changes In 
the nodal parameters, and the mass matrix represents time variation of 
residuals. A solution of (5.1.1) is a possible disturbance or vibration 
mode, the matrix vector X, and the corresponding decay or vibration rate 
suqared, x. These early works also began fueling Interest in the asym- 
metric variety of (5.1.1) which arises when damped vibrations are con- 
sidered. The convention that repeated indices be summed is adopted in 
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Equation (5.1.1) and below. 

The first notable stsp toward an effective strategy for handling 
large asymmetric eigenproblems with the digital computer was by Bauer in 
1957. He devised a scheme in which both dominant eigenvectors - i.e. 
eigenvectors for which the corresponding eigenvalues are largest in 
magnitude - and dominant adjoint eigenvectors - i.e. dominant eigenvec- 
tors of the adjoint eigenproblem 

where X* Is the 1-th component of an adjoint eigenvector and x* Is the 
, corresponding eigenvalue - were constructed simultaneously and the well- 
known property of bl orthogonality (see e.g. Friedman 1956) was exploited. 
Bauer's method was the first to use simultaneous Iteration to find the 
dominant subspace of an asymmetric eigenproblem. The dominant subspace 
Is the subspace spanned by the dominant eigenvectors. However, his 
method has two major shortcomings. The first is that the computational 
work Is doubled by tracking both eigenvectors and adjoint eigenvectors. 
The second is that eigenvectors may be nearly parallel and as a result, 
there can be problems of numerical Instability (Stewart 1976*. Some of 
the practical aspects of Bauer's method and additional Improvements were 
later discussed by Clint and Jennings (1971). 

Work on another computational strategy which would find all the 
eigenvalues and optionally all the eigenvectors of an asymmetric matrix 
had begun somewhat earlier. Von Neumann and Goldstein (1947) first pro- 
posed that the method of Jacobi (1846) should be extended for the reduc- 
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tlon of a general matrix to triangular form. Interest In this Idea did 
not spring up until about ten years later, when Rutlshauser {1955, 1958, 
I960, 1963) developed the LR algorithm, which Is based on repeated factor- 
ization of a matrix Into the product of a lower- triangular (L) matrix and 
upper-triangular (R) matrix. Francis (1961, 1962) Improved Rutlshauser 1 s 
algorithm markedly In the development of the well-known QR algorithm, 
which Is based on repeated factorization Into the product of an orthogo- 
nal (Q) matrix and an upper-triangular (R) matrix (see Stewart 1973 for 
a detailed description). Wilkinson (1965) has called the QR algorithm 
"the most effective of known methods for the solution of general alge- 
braic eigenvalue problems." It remains today the most popular and effec- 
tive method for finding all the eigenvalues of a matrix. Furthermore, 
ft plays a significant role In many computational strategies which seek 
only the dominant sub space. 

With these two pieces In place, I.e. fcuer's and Francis's works, 
the time was right to build up better computational strategies for 
tracking dominant subspaces of asymmetric matrices. In a series of 
papers Stewart (1969, 1975, .1976, 1978) developed a method, referred to 
here as. Stewart's algorithm, for calculating the dominant subspace of a 
real matrix. He drew heavily from Rutlshauser 1 s published algorithm 
(1970) for solving the symmetric elgenproblem, and he Incorporated 
Francis's QR algorithm as well. The (fetalis of Stewart's algorithm are 
explained more fully 1n Appendix A. This algorithm Is used here to 
solve the elgenproblem which arises from the stability analysis. Alter- 
natives are discussed In Section 5.3 and further ideas In Section 5.4. 
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Before proceeding ft is essential to define the stability problem 
to be solved 1n a precise and mathematical way by addressing two issues. 
What class of perturbations are to be allowed? What criterion should be 
used to Judge whether a flow Is stable or unstable to that class of per- 
turbations? The class of allowable perturbations is limited .here by 
four restrictions: {1} they are to be small in amplitude in comparison 
with the base flow; (2) they must satisfy all physical constraints, which 
are local conservation of mass, no slip and no penetration boundary con- 
ditions, and balance of normal stress at a free surface; (3) they must 
be representable in a set of basis functions which are products ''of finite 
element and Fourier bads functions, as described In Section 5.2.; and 
(4) they must disappear far upstream and downstream of the forming zone. 
The first restriction removes the possibility of considering the Interac- 
tions of a perturbation with Itself due to self-convection and to cou- 
pling between free surface elevation and stress components at the free 
surface, but so doing allows great simplification 1n the mathematics. 
In particular, small perturbations have the following properties: 1) 
they grow or decay exponentially in time; 2) they obey the property of 
superposition; and 3) therefore may be decomposed Into Fourier normal 
modes in the transverse direction z, The second restriction Is self- 
evident: perturbations must conserve mass and satisfy essential boundary 
conditions. The third restriction Is a convenient choice of perturba- 
tion representation to make manageable the stability analysis on a digi- 
tal computer. Furthermore,. It is consistent to use the same set of 
finite element basis functions to represent the perturbation as was used 
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to represent the base flow; otherwise the stability analysis .night show 
the base How to be unstable to a mode which the base flow analysis did 
not permit. The fourth restriction is to eliminate from the stability 
analysis consideration of possible modes which originate far upstream or 
downstream of the forming zone. The criterion for Judging stability of 
a flow to the class of disturbances defined above 1s that of asymptotic 
stability, i.e. a base flow 1s considered to be stable If and only If 
all allowable perturbations decay as time Increases, it will be shown 
in Section 5.2 that a baseflow is stable whenever all the eigenvalues of 
a certain elgenproblem have positive real part. 

Sections 5.2 and 5.3 respectively show how to formulate the stabil- 
ity problem In terms of an elgenproblem and how to solve the elgenproblem. 
The Ideas are summarized and compared with some ongoing work of Ruschak 
and of Brown in Section 5.4. A promising variant of Brown's method is 
being Investigated by Basaran and Benner. 

5.2. Formulation of the Stability Problem 

Here the equations governing a small unsteady three-dimensional 
perturbation to a steady two-dimensional base flow are set out and the 
Galerfcln weak form of these equations is constructed. The superposition 
of a small perturbation and a steady two-dimensional base flow can be 
represented by the perturbed velocity and pressure fields: 

u(cl s u(x,y) ♦ eu'lx.y.z.t) (5<2a) 
pie) = p( X ,y) + cp'(x,y,z,t) {5 , 2>2 ) 

The shape of the domain changes because the free surface elevation is 



IS3 

also perturbed: 

TT(e) s h(x) ♦ eh'(x,z,t) . ( 5i2 .3j 

Here overbars Indicate a perturbed field or quantity, primes Indicate 
perturbations to the baseflow, and c 1s a small parameter which Is pro- 
portional to the amplitude cf the perturbation. Secause e is snail, 
only terms that are linear in It need be retained, for higher-order terms 
are negligible by comparison. Formally, linearizing in e is equivalent 
to taking the first derivative of the governing equations and boundary 
conditions with respect toe and then setting e to zero. The strategy 
illustrated In Figure 4.1 to solve the steady base How. i.e. dividing 
the flow domain into three zones, is also used to solve the three- 
dimensional unsteady perturbed How. and thus the equations have to be 
linearized In each of the three zones. 

Figure 5.1 Illustrates the three-zone- strategy for analyzing the 
perturbed flow. Linearization of the governing equations, (1.2.1} and 
(1.2.2), leads to the result: 

3u' 



— - - ne * muu ♦ oi ♦ v . t" (5.2. 4) 



Re ~- = - Re v '(uu' ■•• jj'u) +7 . r 
0 a 7 . „• . 



(5.2.S) 

Here the stress tensor is V « (7u' ♦ 7u' T ) - 7p' because the liquid 1s 
Newtonian and Incompressible and Re s Ud/v is again Reynolds number, u 
Is web speed, d is the upstream gap thickness as Illustrated in Figure 
4.1. and v is thp kinematic viscosity. In the forming zone the base 
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flow velocity, u, Is known as an expansion in finite element basis func- 
tions from the analysis of Chapter 4. In the asymptotic free surface 
and channel zones the base How velocities are solid body translation 
and parabolic, respectively. The acceleration terms In these asymptotic 
zones have been set to zero as an approximation, as Indicated In Figure 
5.1. This approximation Is Justified in the vicinity of a neutral sta- 
bility point when the neutrally stable disturbance does not oscillate 
temporarily. Asymptotic solutions were constructed by perturbation anal- 
ysis for these two flow zones in Chapter 2. These solutions are used 
below to construct boundary conditions at the Inflow and outflow bound- 
aries of the forming zone. 

To expand the velocity and pressure perturbations in three- 
dimensional finite element basis functions In the forming zone would be 
prohibitively expensive for reasonably fine meshes (e.g., a 3 x 13 x 5 
element mesh would Involve about 6500 unknowns). Instead, the device of 
employing Fourier normal mode analysis In the transverse coordinate used 
previously In Chapter 2 is applied here as well, and this results in a 
tenfold decrease In the number of unknowns. Perturbations are repre- 
sented by basis functions which are products of the same two-dimensional 
finite element basis functions that were used to approximate the base 
flow In Chapter 4 and Fourier basis functions. Furthermore, because the 
governing equations (and boundary conditions to appear below) are linear 
In the perturbation, time behavior is exponential. Therefore, a pertur- 
bation of single wavenumber has the following representation: 
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u'(x,y,z,t) » u^'(C,n) • 0(Nz)e" U (s. 2>6 ) 
p'(x.y,i.t) * pj^tt.n) cos(Nz)e" u (5.2.7) 
h'fx.z.t) « hip 1 "!?) cos(Nz)e" U , ( 5 >2 . 8 ) 

where u], pj. and h| are coefficients which determine the fonn of the 
velocity, pressure, and free surface elevation perturbations, respec- 
tively; i, 1 , and p* are the same finite element basis functions 
defined In Chapter 4; the Isoparametric coordinates ($.n) represent the 
spatial coordinates fx,y) as defined by the isoparametric transformation: 

(T(c).y(e)] * ((x 8 .y 0 ) ♦ c(x«^)J f **(«.„) . (5.2.9J 
Here [7 0 (c),y 0 (c)] i are the perturbed coordinates of node 1 which oust 
move with the free surface whenever node 1 lies beneath the free surface; 
M Is the transverse wavenumber and 1s Inversely proportional to the 
transverse wavelength, 2ir/N; D(Nz) 1s a specially constructed tensor 
which serves the same purpose as the vector 0(Nz) defined 1n Equation 
(2.4.9), 

0(Nz) = cos(Nz) (_H + jj) + s1n(Mz) kk ; (5.2.10) 

and A is the exponential decay rate. Summation over the repeated index 
1 is Implied. The exponential decay rate, A, may be complex 1n the 
general case. The real part of A, A R , is the temporal decay or growth 
rate; the Imaginary part, Aj, Is the temporal oscillation rate. When A 
is real the perturbation grows or decays In a purely exponential manner; 
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when ft Is complex the exponential oehavior is nodulated by temporal 
oscillation. The crucial characteristic of A in stability analysis is 
whether its real part is larger than zero, in which case the How is 
stable to the corresponding perturbation, less than zero, in which case 
the flow is unstable to the corresponding perturbation, or equal to zero, 
in which case the How is neutrally stable to the corresponding pertur- 
bation, i.e. the perturbation neither grows nor decays in time. 

Boundary conditions for perturbations in the forming zone are anal- 
ogous to those for the steady base How, as shown in Figure 5.2 (cf. 
Figure 4.3). At the Inflow and outflow boundaries, velocity and stress 
are required to be continuous. At solid surfaces neither slip nor pene- 
tration is allowed, except near the contact line where slip must be 
admitted If the contact line can move. Huh and Scrlven (1971) showed 
that no slip at a contact line is physically inconsistent with movement 
of that contact line. SilHman and Scrlven (1978) showed that allowing 
sane slip near a contact line, either by prescribing a slip coefficient 
or slip velocity profile, removes the stress singularity which would 
otherwise result. When a contact line can be perturbed, i.e. when it is 
not pinned at a geometric or compositional discontinuity, it is necessary 
to permit some slip. This is done here by prescribing a slip velocity 
profile that Is convenient in the finite element formulation. Boundary 
conditions at the free surface are the following: (1) the kinematic 
boundary condition states the free surface must deform at a rate consis- 
tent with the arrival or departure of liquid there; (2) because the gas 
phase is considered Invlsdd as In Chapter 4, the shear stress boundary 
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condition states that a perturbation of velocity and free surface eleva- 
tion must oaintaln zero shear at the Interface; and (3) the normal stress 
boundary condition states that change In the nonial viscous stress and 
pressure at the interface must be balanced by a change In the capillary 
force due to the mean cv-vature. H. Additional boundary conditions at 
the perimeter of the fr« -.urface are that the interface be continuous 
at the outflow boundary, both in elevation and slope, and that either 
the contact line or the contact angle not be disturbed. The boundary 
conditions at the contact line correspond respectively to cases where 
the contact line is pinned or free Jo. slide along the solid boundary. 
The free contact line case is developed here; the pinned contact line 
case is Identical except that the free surface elevation at the contact 
line is an essential boundary condition, and thus the contact line term 
in Equation (5.2.33) below can be omitted. Boundary conditions at the 
ed S e surfaces, I.e. planes of constant r, need not be written explicitly 
because the domain is chosen to be one wavelength, 2»/N, In the trans- 
verse dimension so that periodicity conditions apply. The boundary 
terms at the two edge boundaries cancel each other as shown below. The 
set of boundary conditions described above is: 



V -4 ,a0 (5.2.1J.) 

i-<T-S>-£ (5.2.12) 

h i h . . Outflow boundary 

h ' "V " 0 (5.2.13) 

d(h' - hi)/dx » 0 (5>2a4) 
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i* ( i' ~& 30 (S.2.1S) 



I" (5.2.16) 



u 1 xn a 0 



Inflow boundary 
Solid surfaces 



(5.2.17) 
(5.2.18) 



u'x„ 3 -4 L Q( £ ,jc Solid (SAW) 

(5.2.20) 



Solid surfaces 
u 1 • n » 0 rtear contact 1 Ins 



(5.2.21) 
(5.2.22) 



ah'/3t * u' • v(y-h) - u • vh' 
0 * [(7 • ^T) x n]' Free surface 

0 - [Cn - T) - n]< . c^K' (5 ^ 23) 

e * B 0 'I Contact line (5.2.24) 

The unit normal n at the boundary 1s Illustrated In Figure 5.2. u 1 Is 

CL 

the slip velocity at the contact line and Q(e) is 4 quadratic function, 
of the Isoparametric coordinate 5, as explained tela*. c« = v u/« is 
again the capillary number, u 1s* dynamic viscosity, U Is web speed, a 
Is surface tension, and h" Is the change In the mean curvature due to a 
perturbation of the meniscus elevation, h'. This set of boundary con- 
ditions. Equations (5.2.11)-(5.2.24), applies generally to small pertur- 
bations to coating Hows. However, the physics governing the inmediate 
vicinity of a wetting line. I.e. the contact angle and the flow nearby, 
1s not yet firmly established and thus Equations (5.2.19) and (5.2.24) 
may be modified or replaced as understanding Increases. 
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The Galerkln weak forms of the unsteady governing equations were 
derived In Section 4.2. Galea's aethod requires that the basis func- 
tions which appear In Equations (5.2.6MS.2.8) be chosen as weighting 
functions, and thus the wight**, residual equations for the perturbed 
flow field are: 

■ I R6 "tr 1 ' I* 1 "*) Win)) ^llff fcdndz 
3 / {Re 7 . ©eJuUJ].[^{c,n) 0(Nz)] 

♦ T(c):» n) 01Hz)] ill^iii d5dnd2 

~ _ 3(C,n) 

- / JU) • T(e) • [#J{ €t ,) o(Nz)l il£l dy dz 
(«0-*Q J " 3 Y 3z 



• - / n(e> . T(e) . [*J( 5l „) 0(Nz) iHM^cli 

6fl. " - »U.nl atonaz 

0668 (5.2.25) 

o - / v • u (e ,*J (c . n , »P»;{goi d5dnd2 f (s>2<26) 

where 0 represents the three-dimensional donaln. 6a represents the bound- 
ary of o and has five distinct parts: 

50 3 M SS + ««FS + 50 IM * «0UT + «EBBES ' (5.2.27) 
These respective parts are solid surfaces, free surface, Inflow boundary, 
outflow boundary, and edge boundaries, i.e. planes of constant trans- 
verse coordinate z as illustrated In Figure 5.2; Kc) is the perturbed 
arc-length along the boundary at a plane of constant zfc(e) i s the per- 
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turbed arc-length along the boundary at a plane of constant streanwise 
coordinate x; and a[x(e)J(e J]/3(5.n) fs the Jacobian of the Isopara- 
metric map. As mentioned above, a convenient choice of domain width Is 
az » 2b /H. one wavelength, because- the boundary Integral over the edge 
boundaries In Equation (5.2.25) drops out by cancellation of the terms 
at opposing edges. 

The next step in the derivation Is to linearize the governing equa- 
tions. (5.2.25) and (5.2.26), In the parameter*. This Is a complicated 
step because the nodal parameters, u,. p r and h,. and the nodal coor- 
dinates (x 0 ,y 0 ) i all depend on e. The nodal coordinates do so because 
as the free surface deforms, the nesh lying beneath also deforms. This 
affects the acceleration tern In Equation (5.2.21) because points of 
constant e and n are convected through space and therefore (cf. Kheshgl 



1982) 

»"<e) | „ jufejl , ahj(e) ax(c) I 

(5.2.28) 

The last tern In this equation Is evaluated as follows: (1) the first 
factor Is the rate of change of the free surface elevation parameters; 
(2) the second factor, a vector, 1s the rate of change of the nodal 
coordinates with the free surface elevation parameters, which depends 
entirely on the mesh construction; (3) the third factor, a tensor, Is 
the rate of change of global coordinates. (x,y). with nodal coordinates 
at fixed isoparametric coordinates, !c,r.». as determined by Equation 
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(5.2.9); and (4) the last factor, a tensor, is the gradient of the veloc- 
ity vector. Linearizing and solving for au'/atl i M ds to 



au' ^ au/ 
3t U.y »t 



ax! , . 
♦ xcosOizlhi -3. . JUL . v u /c , , 0 » 
«' n J9h j %l(l.n) " (5 • 2 • ^9, 

The gradient terns in Equations (5.2.25) and (5.2.26) also depend on 
both nodal parameters and nodal coordinates. This is more easily seen 
from the equation for the gradient in terms of the isoparametric coor- 
dinates, viz. 

3x - 3y - 3z- 

, c a[7( e ),7( e )] r i ff ay( e ) j_ s^e) a 
»(«.n) J ^ an a? ' aT aT^ 

+ f-Ulii i_ + i£f£i 3 m «, a . 

3n a? 3g 3V iJ il^- (5.2.30) 

The last equality results from the inverse of the isoparametric map. 
Equation (5.2.9) is needed to evaluate the derivatives of global coordi- 
nates, [xtc),y{c)], with respect to the isoparametric coordinates, {?,„). 
The linearized form of (5.2.30) is fairly complicated, and fortunately 
need not be computed because the inverse Jacobian of the isoparametric 
««P in (5.2.30) cancels the Jacobian in (5.2.25) and (5.2.26). The 
linearized form of the Jacobian times the gradient is (Safto and Scriven 
1981): 
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u an ac 3? 3n Ji 1 3n a« aT j^ JiJ ' (5 ^ 3l) 

Equations (5.2.25) and (5.2.26) can be Unearllzed with the aid of Equa- 
tions (5.2.29) and (5.2.31). The results are exemplified by the deriva- 
tive of the J-th x-mcroentum equation with respect to the 1-th x-veloclty. 
unknown, uj, when the element that contains uj borders on the free sur- 
face; 

0 = - Re X / *S J cos 2 (Nz)[3(x t y)/3U,n)]dCdndz 

♦ / (Re [4> j 3u/3x + u3* f /3x + v3$ f /3y]A J 
Q 

+ [(23<frVax)(34> J /3x) + o+'/ayHa+'/ay) 

- M 2 *V]} cos 2 (Mz) l f 3{x f y)/3(€ l n)]dCdn(lz 

- / (n*[2l3« f /3x + J 3^/3y] cos 2 (Nz)(3s/3C)<«dz (5.2.32) 

The boundary conditions are administered in the same way as In 
Chapter 4. The adherence boundary conditions In (5.2.17)-(5.2.20) take 
the place of the momentum equations (5.2.21) weighted by basis functions 
associated with boundary nodes, and thus the surface Integral over on., 
need not be evaluated because each trial function associated with all 
other than solid boundary nodes 1s Identically zero at the solid bound- 
aries. Equations (5.2.17) and (5.2.18) are used everywhere along the 
solid surfaces except In the element adjoining the contact line, where 
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(5.2.19) and (5.2.20) are prescribed ind Q(e) is chosen so that ft Is 
unity at the contact line and it and its first derivative disappear at 
the .ipstreasi edge of the element. The shear stress boundary condition, 
Equation (5.2.22), is imposed as a natural boundary condition in Equation 
(5.2.255 by setting 

n.T-nnn:T (5.2.33) 

in the free surface boundary integral. The y component of each momentum 
residual, i dotted into (5.2.25). 1 S replaced at free surface nodes by 
the kinematic equation. (5.2.21). The linearized normal stress equation 
Is enforced along the free surface: 

0 = / {[2n'n:T* n n:T']cos(Hz)p^c)<a"^^?[co$(Hz)p^Oj^*dz 
FS 3 5 

♦ / [n n:T cosfNzJp 1 ^) + Ca'VvfcosfNzJp^}]} iil dedz 

«0 FS - - 3 € - 

" U '\t k''!L c «(Nz)p f (5)]ds , (5.2.34) 

««n FS 

The last tern is the line Integral around the perimeter of the free sur- 
face boundary and m is the outward pointing normal to the base plane, 
i.e. xz-plane. The line integrals along the two edges of constant trans- 
verse coordinate i cancel each other and thus the last terra simplifies to 



" Ca _1 j5 [n'.tn cosfMzJp^Ojds =» - Ca' 1 * [n'-i cosfNzJp'ujtfhr dz 
«6B_ S 64a ^- - 'OUT 

where CI and OUT refer to the contact line and outflow boundary, respec- 
tively; the normal to the free surface is 
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n s n ♦ e n' ♦ ... 



, ' ' h ^-lVl)cos (Hz)^(Uh^h'3ln(N2)lc 

" / T7T ? € ^n?P72 = • (5.2.36) 

X x 

Boundary conditions (5.2.13) and (5.2.14) combine to fonn a Robin con- 
dition at the outflow boundary: 



n' . 1] 0UT * . „{ N )h'/[l ♦ o 2 (0)(h - hj 2 ] 3/2 



(5.2.37) 



Here «(0) and «(N) are the slowest decay rates for perturbations in the 
asymptotic free surface zone when wavenumber Is 0 and N, respectively, 
as calculated In Chapter 2. Equation (5.2.38) 1s substituted In Equation 
(5.2.35). Equation (5.2.24) with the aid of the trigonometric relation- 
ship 

e « 180- . tan' 1 h x - tan' 1 hj . (s . 2 . 38) 

is also imposed as a natural boundary condition in Equation (5.2.35). 
Here h s Is the elevation of the solid die above the web. Boundary con- 
ditions (5.2.15) and (5.2.16) combine to give a Robin condition at the 
inflow boundary and (5.2.11) and (5.2.12) combine to give a Robin con- 
dition at the outflow boundary. These are Imposed as natural boundary 
conditions in the linearized form of (5.2.25): 



/ !• r . 0(Nz)* 1 (e.n)&*Hfc 
40 IH " " 

■/ n 8. «■ • 0(Mz)^( C .n)£dndz , (5 . 2 . 39) 
IN 
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- / i- T' • OfNzH'u.n) ^dndz 
Mu.- - - dn 




'out ~ 



(5.2.40) 




(5.2.41) 



As {(i(2|f -p')/u' + i^ u . ♦jciul/u'Ji 



(5.2.42) 



Subscripts C and F Indicate that the quantities Inside the brackets per- 
tain respectively tn the upstream channel and downstream free surface 
zones. Chapter 4 demonstrates that a Robin condition is better than 
Neumann or Dlrlchlet conditions at an open-flow boundary (especially at 
an outflow boundary) for the two-dimensional flows considered there. 
The Robin conditions in Equations (5.2. 39}-(5. 2.42) are three-dimensional 
analogues of the one in Equations (4.2.19) and (4,2.22), and so they 
extend from two-dimensional flows to tfcrss-diaensional ones the idea of 
matching asymptotic zones to forming zone. 

The final step In the derivation of the equations governing a per- 
turbation to a steady two-dimensional flow Is to integrate over the 
transverse coordinate, z. This turns out to oe trivial because every 
term in the governing equations, (5.2.25) and (5.2.26), and in the normal 
stress equation, (5.2.34), is weighted by either cos 2 (Mz) or sin 2 (Nz) 
(cf. Equation (5.2.32)) and these both Integrate to the same constant. 



168 

The final form of the resulting equation set Is an asymmetric generalized 
eigenproblem, 

J ki ^i 3 X \i «1 ■ (5.2.43) 

as exemplified by Equation (5.2.32) and Illustrated by Figure 5.3. The 
symbols In Figure 5.3 are defined as follows: 

Mx 

i Right side of 1-th x-raanentum equation, 1_ • (5.2.25) 

Mv 

Rj J = Right side of 1-th y-momentum equation. j_ • ;5.<:.25) 
Mz 

Rj = Right side of 1-th z-racmentum equation, £ • (5.2.25) 
K 

R t s Right side of the 1-th kinematic equation, (5.2.21) 

rJ h 1-th continuity equation, (5.2.26) 
NS 

R,| h 1-th normal stress boundary condition, (5.2.32) 

M 1j 5 ^" th element ° f Inertial mass matrix, i.e. the derivative of 
the 1-th left side of the x-momentum equation, 1 • (5,2.25). 
with respect to u % . ~ 

M^j s ij-th element of x-convective mass matrix, 1.e. the derivative 
of the f-th left side of the x-moracntum equation, 1 • (5.2.25). 
with respect to hi - * 

Cy 

M 1j ' ^" th e1ement of y-convective mass matrix, I.e. the derivative 
of the 1-th left side of the y-moraentum equation, J • (5.2.25) 
with respect to h j - * 

H^j 5 IJ-th element of kinematic mass matrix. I.e. the derivative 
of the 1-th left side of the kinematic equation, (5.2.21). 
with respect tu hj 



The adherence boundary conditions, (5.2.17) and (5.2.18), result In the 
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Kronecker delta entries in Figure S.3. The eigenvector In Equation 
(5.2.43), the i-th element of which 1s denoted by q { , represents the 
velocity, pressure, and free surface unknowns, I.e. 

W 5 fu l V v l V w l V p l fra' "i.-mHJ 1 . 

Here n Is the number of velocity nodes, m Is the number of pressure 
nodes, and i the number of free surface elevation unknowns. Each elgen- 
solutlon Is a possible perturbation nude and Its decay rate. A base 
flow is therefore stable If each possible mode has a positive decay rate 
at every value of wavenumber, Jl. Means to solve this problem are the 
topic of Section 5.3. 

Ruschak (1980) and nore recently Klstler and Scrlven (1982) have 
described an alternative method for introducing the boundary conditions 
at the free surface. They suggest that both normal and shear stress 
conditions should be used as natural boundary conditions for the weak 
form of the momentum equations. In their formulation, the free surface 
boundary Integral In Equation (5.2.25) would be replaced by 

■ Ca" 1 / 2H n • 0 cos(Nz» f (c,n) S to . (5<2#44) 
6Q FS " H 32 

The surface divergence theorem (Weatterburn 1927) can be applied to the 
right side to get: 
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/ 2Hn. D CostNz}*^,,,) ilde ^dz 
40 FS H n 

...1 



\L ? f - (MZ) * U,n,dS • (5.2.45) 

Here 7 $ is the perturbed surface gradient and m Is the perturbed outward 
pointing tangsst on the perimeter of Uic free Surface. Evaluation of the 
last terra in this equation is analogous to evaluation of the perimeter 
integral of the normal stress condition in (5.2.35). The surface diver- 
gence terra is evaluated by using the definition of the surface divergence 
(Weatherburn 1927). 

5.3. Methods of Eiqena nalysis for the Stability Problem 

Tha festures of the elgenproblem defined in the previous section 
are the following: (1) the matrices are large in size, about 650 by 650 
for the raesh used in Chapter 6, but are sparse and banded; (?) both the 
Jacobian and raass matrices are asymmetric; (3) the mass matrix is singu- 
lar because continuity, solid wall velocity, and normal stress constraints 
on the perturbation contain no time derivatives (see Figure 5.3). These 
features are crucial to the choice of an algorithm to solve the eigen- 
proble- The first feature rules out as too costly any algorithm which 
coop. : ; the eigenvalues or requires full matrix storage. The second 
precipes most of the available algorithms because they cannot handle 
the asymmetric elgenproblem. The third feature, that the mass matrix is 
singular, eliminates the possibility of Inverting that matrix. Thus it 
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is not possible to reduce the generalized elgenproblem to a simple eigen- 
problem of the form: 

Because of the above restrictions there are only a few viable methods 
to perform the eigenanalysls. Three possibilities are worthy of consid- 
eration: (1) Mewton's method; (2) Inverse iteration; and (3) inverse 
subspace iteration. The pros and cons of these methods are described in 
the following paragraphs. 

Newton's method is the method of choice for solving the nonlinear 
problem of Chapter 4. i.e. solving the two-dimensional steady baseflow, 
and therefore it would seem te be a worthy candidate here as well. To 
apply Newton's method. Equation (5.2.41) Is rewritten as a z±z of resid- 
ual equations: 

Rj » (J^ - W^Jq, , 1 - I.....L, (5.3.1) 

where L is the number of degree's of freedom In the basis set. An addi- 
tional constraint Is needed to normalize the components of the eigenvec- 
tor, q^. A convenient choice Is: 

\*1 "J (1 ' *V (5.3.2) 

Applied to Equations (5.3.1) and (5.3.2), Newton's method gives the 
following Iteration formula: 
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. fl (k) 
q 1 


0 



r < m) 



Jk) 



Jk+D . x (k) 



« l J i 



„(k) 
K t+1 



k ■ 1,2,..., 



(5.3.3) 



where k is the Iteration count and q{ 0J and x<°» is an Initial guess. 
Newton's method has the advantage of quadratic convergence of tee update 
vector, (6q/ k) . ix M ), to zero with iteration count k. However, the 
drawbacks to using Newton's method to perfonn the eigenanalysis are many 
and serious: (1) the linear system of equations in (5.3.3) would have 
to be solved several times to find each elgenpair; (2) complex arithmetic 
would have to be used extensively if a complex elgenpair were sought, and 
this would double the memory requirement for Jacobian and eigenvectors 
and thereby increase computation costs significantly; and (3) the elgen- 
pair found at convergence of Newton's method would depend strongly on 
the Initial guess. The last point is the most severe drawback because 
the goal of the analysis fs to find the most dangerous modes, 1.e. the 
ones wh-. >cay most slowly or grow most rapidly in time. With Mason's 
method the- is no assurance that the elgenpair found 1s the most danger- 
ous unless a very good Initial guess Is known. 

The second candidate, Inverse Iteration, is simpler to implement 
than Newton's method. At each stage of the iteration the linear set of 
equations defined by 
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(5.3.4) 
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where qj 1 ' Is a.jain an initial guess and qj k+1) 1s the (k+l)-th Iterate 
prl;r to normalization, i.e. 

q{ k+1) =qj k+11 /. qj k+1) . (S.3.5) 

where i • J Is a convenient norm. The corresponding eigenvalue 1s given 
by 

This algorithm has been shown to converge linearly to the elgenpalr which 
has smallest absolute value, provided one exists (Stewart 1973). The 
chief advantage of inverse Iteration over Newton's method is that inverse 
iteration 1s more likely to converge to the dominant mode, but on the 
other hand 1t may converge slowly or not at all depending on the spectrum 
of the problem, I.e. the distribution of eigenvalues (see Appendix A for 
further discussion). An additional advantage of Inverse iteration (also 
Inverse subspace iteration) over Newton's method 1s that Hood's method 
(1976), which is used to solve the base How equations as described in 
Chapter 4, can be used to solve (5.3.4). 

The third candidate, simultaneous Inverse iteration, 1s Hke simple 
inverse Iteration except that several eigenvectors are tracked simulta- 
neously. Algorithms to perform simultaneous Inverse Iteration are of 
necessity more complicated than algorithms for simple inverse iteration, 
but they offer Improvements. Some of these are: . (1) the convergence 
rate can be markedly increased (see Appendix A for details); (2} complex 
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eigenpairs can be found easily without need for explicit complex arith- 
metic; (35 dominant modes can be found simultaneously instead of sequen- 
tially as in the simple inverse iteration scheme. The algorithm of 
Stewart (197G. 1978) is the most elegant of the alternative simul taneous 
Inverse Iteration algorithms for asymmetric matrices (see Bauer 1957. 
Clint and Jennings 1971. and Dong 1977). Appendix A discusses at length 
Stewart's algorithm and Its Implementation. 

Inverse subspace Iteration aay not always find the »st dangerous 
modes, however, it finds the smallest set of eigenvalues In absolute 
value. i.e. the ones for which | A | * /a* + X J 1s ^ miUe ^ m ^ 
most dangerous modes are the ones for which A R is least, i.e. A R is most 
negative. At some safe parameter values the flow should be stable and 
thus each A R should be greater than zero. When marginal stability Is 
approached, for example by Increasing capillary number (see Chapter 6). 
one or more of the A R becomes nearly zero. If a r passes through zero ' 
and becomes negative, then the solution bifurcates and loses stability. 
If the eigenvalue passing through zero is real. 1.e. ^ = 0 . then In- 
verse subspace iteration cannot miss the bifurcation point because |x| - 
0 there. But. the bifurcation point would be missed if Aj were large in 
comparfson with the magnitudes of other eigenvalues, I.e. if the domi- 
nant mode had , high temporal frequency. Computation of the onset of 
high frequency modes of Instability requires a special trick from the 
methods of conformal mapping. The bilinear transformation is known to 
transform generalized circles in one complex plane into generalized 
circles in another complex plane (Churchill et al. 1974). By choosing 
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the parameters of the capping appropriately, a lino can be transformed 
into a circle, and the half-plane to the left of the line into the 
interior of the circle, as illustrated in Figure 5.4. A bilinear trans- 
formation that does this Is 

U3 2?T ' (5.3.7) 

where d 1s the x-lntercept of the vertical line shown in Figure 5.4. 
The Inverse transformation, also bilinear, is 

X ' (5.3.8) 

Replacing the eigenvalue X In Equation (5.2.40) by (5.3.8) and rearrang- 
ing gives 

J jl1i a * (2d H J1 * Jj f >q f . (5.3.9) 

which fs another generalized elgenproblem, but with » as eigenvalue. By 
choosing d carefully, the most dangerous modes can be found by solving 
(5.3.9), However, the choice of d is crucial; It oust be so chosen that 
the number of eigenvalues in the half-plane to the left of d Is exactly 
the same number as those sought in the eigenanalysis. If more are pres- 
ent, the ones found by the eigenanalysis would not necessarily be the 
scst dangerous; if fewer are present, Inverse subspace Iteration would 
not be likely to converge. Convergence would be lost because of the 
structure of the spectrum In the u-plane. The distance of each eigenval- 
ue from the origin in the y-plane is found from Equation (5.3.7) to be 
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M-IM/Z^r.- 2d ' 2 + x I . (5.3.10) 

Eigenvalues at Infinity In the A-plane are mapped Into ones which are zt 
a distance of unity from the origin 1n the y-plane. Unfortunately, cany 
such eigenvalues exist, one for each continuity, adherence, and normal 
stress constraint on the perturbation. As explained in Appendix A, 
Inverse subspace iteration does not converge when the dominant subspace 
is not uniquely determined* Furthermore, even if d were chosen properly, 
inverse subspace Iteration might converge more slowly for the mapped 
eigenproblem in (5.3.9) than for the original eigenproblem in (5.2.43) 
because large eigenvalues 1n the X -pi arte would He near the unit circle 
In the u-plane, as seen from Equation (5.3.10). (The effect of the spec- 
trum on the convergence rate of Inverse subspace Iteration 1s quantified 
in Appendix A.) The bilinear mapping, then, should be used only with 
caution and when temporally oscillatory modes of instability are sought. 

8ecause Inverse subspace Iteration Is an Iterative technique for 
solving a nonlinear problem, an Initial guess is needed. The best w^y 
to choose an initial guess 1s to use Information from nearby solutions, 
provided one or more are available. If tne dominant subspace has been 
calculated at an. old set of parameters, that subspace can be chosen as 
an l.iltlal guess of the dominant subspace at new parameter values, a 
technique known as zeroth-order continuation. A better Initial guess 
can be constructed by using first-order continuation. In which the rate 
of change of the dominant subspuce with respect to the parameter or 
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parameters to be changed Is calculated at the old point and used to con- 
struct an extrapolation in Euler-Hke fashion to a new point. However, 
first-order continuation is more difficult and more costly to implement 
than zeroth-order continuation, and 1s therefore preferred only wnen the 
savings gained (n the iteration process because of improvement of the 
initU guess outweighs the additional cost needed to construct the ini- 
tial guess. In the case of the generalized asymmetric elgenproblem, the 
advantage of first-order continuation is doubtful. For example, to cal- 
culate the rate of change of one eigenvector to a parameter, say s , the 
derivative of Equation (5.2.41) with respect to s is calculated: 

( 'V as - V«1 * CJ Jf - XMj^yas . o . (5,3.11) 

The derivative of the Jacobian matrix with respect to the parameter, 
3J j1 /3s ' has to be calculated and the Equation (5.3.11) used to find 
ayas and a/3s. If the elgenproblem were symmetric, 3x/3s could be 
found first by talcing the inner product of Equation (5.3.11) with the 
eigenvector and taking advantage of the fact that 

(J j1 - X V ^V 35 ^ B < J ji ' *« Jf >q, Oqj/ss) = 0 , (5.3.12) 
so that 

» OV 88 ^ Wl"j- ' (5.3.13) 

When the eigenproblem is asymmetric this simplification could be made 
only if the adjoint eigenvector were found, an unecononXcsl option. The 
other possibility is to solve (5.3.11) plus an additional constraint, say 
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qfOqf/as) a 0 , (5.3.14) 

simultaneously for 3qj/3s and 3X/3s. This also would be difficult 
because the system of equations to be solved would not be banded and 
therefore subroutine FRONT (Hood 1976) could not be used to find 3q^/3s 
and ax/ss* Zeroth-order continuation Is therefore preferred to generate 
an Initial guess, or when several steps In the same parameter are taken, 
a secant extrapolation could be used. 

5.4. Discussion 

A general method for finding the stability of a steady, two-dimen- 
sional, viscous free surface flow has been described. Furthermore, the 
dominant eigenvectors represent the most dangerous disturbance modes, and 
thus both the onset and mode of Instability can be predicted. Although 
the derivation of the stability equations has been aimed toward coating 
flows In particular, the methods and strategies have obvious general iza- 
tlons to practically ar*y steady two-dimensional viscous flow, with or 
without a free surface. Allowable perturbations are two-dimensional 
(M=0) or three-dimensional, the basic restriction being that they must 
be representable in basis functions which are products of finite element 
and Fourier bases. The methods for solving the resulting eigenproblem 
are general in the sense that ar\y mode of disturbance, even if 1t Is 
time oscillatory* can be tracked. 

Ruschak (1982) has recently developed an approximate method that 
should be useful for creeping flows, I.e. for flows when Reynolds number 



181 



is sufficiently small or for predicting the onset of Instablli^ to;~. - 
disturbances that do not oscillate in tim.,* assumes that acceleration 
is negligible in ^ nonentam equations, i.e. ne ^ates the quasi-steady 
assumption; an ^su^tio^ that Is entirelyjustified only when Reynolds . 
number is sufficiently small. This assumption allo-s him; to elimi.nate - 
each of tte velocity and pressure unknowns from.the eigenproblem. and 
thus to-sta tically con*nse" the ^ original eigenproblem to a much .smaller : 
one. The small elgenproblem can then be solved efficiently by . the qz c 
algorithm (Holer and Stewart 1973). for which IMSL and EISPACK subrou- 
tines are available, ^chafe's approximation is n, 1te attract! veW L 
low Reynolds number flows, but it is not reliable at even moderate 
Reynolds number, in contrast to the one used here.' " 

Brown (1982, has investigated the loss of stability of stagnant 
liquid between two discs heated from below, a Rayleigh^effreys' probl«o 
(sometimes mistakenly called a Benard problem) by finite element analy- 
sis. He considered only axlsymnetric disturbances, however, and thus 
needed no more than a two-dimensional representation. It (s well known 
that the stagnancies layer goes unstable to disturbance modes In 
which convective cells form as adverse temperature gradient is increased 
(Chandrasekhar 1961). Furthermore. Chandrasekhar (1961) proved that 
the convective cells at onset are steady rather than time periodic 
because the decay or growth rate of each disturbance mode is real, i.e. 
the eigenvalues are all real. Finite element stability analysis of the 
Raylelgh-Jeffreys problem leads to an eigenproblem of the same form as 
(5.2.43). Because the Jacoblan ™atHx in this oigenproblem has only 
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real eigenvalues. Brown is able to detect the onset of an Instability by 
examining the sign of the determinant of his Jacoblan matrix. By doing 
sc he avoids the unpleasant task of solving the elgenproblem. However, 
his method could fall If applied to problems for which the eigenvalues 
are not known. a prion to be real. Furthermore, Brown's method might 
have difficulty m resolving closely spaced eigenvalues as they cross 
the Imaginary axis. 

It Is unlikely that the ultimate methods for computer-based stabil- 
ity analysis have yet been found. The work presented here is the first 
of Its kind, although It Is paralleled by those of Ruschak and Brown. 
Furthermore, good computational strategies for handling asymmetric elgen- 
problems have only recently been devised, and even these are still 1n 
their infancy. Giant steps in Improvement of numerical and analytic 
techniques are almost certain to be made In the future, and will enable 
coating flow stability analysis to be done with much more ease and 
simplicity than at present. 
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CHAPTER 6 — STA BILITY OF SLOT Alffl KNIFE COATIMfl .■ «...v^: 
6.1. Introduction ^' "''* 

. . ^ «^»ds ^scH'bed in Chapter 5 were designed to analyze ihe^ta- 
bflity of ftMv? S^iSf«C viscocapil^ry nois..' and furthemore 
to understand the most dangerous 1 disturbance modes, i.e. the patterns, of 
flow that decay.»ost slowly or grow > most rapidly. In this chapter; the. . 
power of these methods is illustrated by analyzing the stability of slot 
fl.vJ knife coating! :-, :'5: . " 

Steady two-dlmenslonaf slot and knife coating flows were treated ■ 
in Chapter 4 by the Salerkin/finite element rcthod. As discussed in 
Chapters 1 and 3, it i S ^^Icnown that these flows can be unstable to . 
three-dimensional disturbances that are transversely periodic, an instal 
b,11ty that ^VC«»only known « ribbing . The parameters that determine 
the stability of the base flow are capillary number; Ca 5 ,u/ a ; Reynolds 
number, Re , Ud/v; dimenslonless flaw rate. Q ; h./d; die shape as deter- 
mined by the exponential factor 1/B defined in Section 4.2; and contact 
line physics, i.e. whether the contact line is pinned at a geometric or 
compositional discontinuity or free to move along the die with some pre- 
scribed contact angle e. The operating parameters and liquid properties 
are as follow: „ i s the dynamic viscosity of the liquid; U,is the web 
speed; a is the surface tension of the liquid/gas Interface; h_ is the 
final uniform film thickness; and d is the upstream gap width oetween 
the stationary die and moving web. The effect ef gravity is considered 
to be negligible in comparison with the effects of viscous stress. 
Particulars of the method of stability analysis are described in 
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Section 6.2. The Influence of liquid and operating parameters on the 
stability of knife coating Is Investigated In Section 6.3. In Section 
6.4 the results are summarized and compared with the approximate rasults 
based on lubrication theory which Mere described In Chapter 3. 

6.2. Stability Analysis 

This section describes the particulars of applying the methods of 
Chapter 5 to the stability of steady two-dimensional slot and knife 
coating flows similar to those investigated In Chapter 4. Attention is 
concentrated on the issues which were not considered in Chapter 5, namely 
the subdivision of the domain into elements, startup and continuation 
strategy, and approximate computation costs. 

The domain is tessellated into elements in exactly the same way as 
in Chapter 4, i.e. by erecting spines orthogonal to the moving web along 
which nodes are placed according to preasslgned weighting factors. The 
node locations define element shapes through the Isoparametric mappir" 
set out in Equation (4.2.4) and illustrated in Figure 4.3. Thus, nodes 
lying on or beneath the free surface adjust in proportion to the local 
free surface eleva ;io«. Tarttizrzsrz, the node at the three-phase con- 
:act line Is free to move along the solid die when the contact line is 
not pinned. Thereby, the spine which Intersects the contact line is 
required to move upstream or downstream with the contact line. The two 
spines immediately upstream and downstream of the contact line are 
required to maintain a position midway between their nearest neighbors. 
Because only small amplitude disturbances are considered, the contact 
line moves only infinitesimally so that no overtaking of spines by other 
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sp1fle $ ;is ;? o^b^..; ^ spt^ were allowed; * cross,, elements -which* - 
have wsj^ve, areas would be. created!. • . ^i^ : >v^- . - . 

danlnant subspace, and^then us! W Stewart'silnye^: subspace 1 Iteration^ 
algorithm (described In Appendix A) until the jnaxiraun error in the e(ge n 
problem residual j was Jess than 10" 8 . o„ce an elgensoScC of ita ^ 
startup ««^;3e*rate^^ 

described In Chapter Sf 'to tract tte to nart : subspace" as**; paraneter 
was varied. Flrs^ a point of neutral stabl 11^, u e . a poi „t at which' 
the dominant decay; rate W 2 erof was sought^ using secant^extr^laV'' 
tion to vary; capillary ^rjai o&rparaSers SSxJdf jhVnt 
the parameters «« varied one-by-o«e> explore their effect x>n stabll^ 
ity. Whenever a new set'of parameters was Investigated, the following . 
procedure ws used: ft) a base flowW calculated;'^' tte^coblan and 
mass matrices were computed fo^ a specified wavenumber and Hood's ^(1976} 
subroutine FROHT was used to execute the forward elimination" part' of 
Gaussian elimination; <3) Stewart's algorithm In conjunction with Hood' s 
(1976) subr6ut1.«^S0L w«% e d to'flnd the dominant subspace for the ? 
given wavenumber; and (4) steps (2) and (3) were repeated for several 
different XrSj^^y^ ; dominant ^space was calculated . 
for five wavenumbers saT that growth rate curves could'oe fared 1n by eye 
(with a French curve). ' " 

The cost of calculating one stability point. I.e.' finding the doml-, 
nant subspace at a given >t of "parameters, ranged from^about one to " 
".ore than ten tlnies the cost of a base now calculation. I.e^ from about 
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55 to 1CCC cpu seconds on the University of Minnesota Cyber 730 (equiva- 
lent to 35 to 600 seconds on the University of Minnesota Cyber 74). The 
cost was relatively snail whenever the ratio oi the magnitude of the 
second eigenvalue to the magnitude of the first eigenvalue was of order 
1C or larger; It was relatively large when the same ratio was of order 
unity, and especially at start-up when no previously calculated dominant 
subspace was available for continuation. To cut down on computation 
costs, It was therefore necessary to use a continuation strategy and to 
follow neutral stability curves. 

6.3. Stability of Slot and Knife Coating 

Figure 4.12 shows a set of streamline contours for a knife coating 
flow at values of ths opcrstlng parameters which are typical of indus- 
trial coating devices and for the case where the contact line Is free to 
slide along the solid boundary. 

Some streamline contours of the most dangerous two-dimensional mode 
of disturbance to the flow Illustrated In Figure 4.12 are shown in Figure 
6.1. This disturbance decays In tlsa, and thus the base now 1s stable 
to It. The streamlines have an eddy-Hlce appearance, except that they 
terminate at the free surface. {Though one end of each streamline 
appears to Intersect the upper solid die. this Is an error in the plot- 
ting routine caused by insufficient resolution. In reality, they inter- 
sect the free surface close to the contact line.) The Intersection of 
the streamlines with the free surface indicate flow normal to the Inter- 
face., and therefore the interface location must be changing with time. 
The recirculatory pattern accompanies displacement of the contact line. 
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For example, if the contact line 1s moving tangent to the die In the up- 
stream direction, liquid upstream of the contact line must be displaced. 
Because the upstream flow rate i< specified, the displaced liquid in 
order to conserve mass, must move downstream of the contact line and 
thereby increase the How rate Just downstream of the contact line. The 
result is that the meniscus must bulge outward. 

Figure 6.2 shows how the decay rate of the most dangerous three- 
dimensional disturbance varies with the wavenumber of the disturbance, 
I.e. with 2tt divided by its dlmensionless wavelength. The two-dimen- 
sional disturbance portrayed in Figure 6.1 corresponds to the limit as 
wavenumber approaches isro, I.e. as transverse wavelength increases 
without bound. Disturbances become increasingly three-dimensional as 
wavelength decreases and wavenumber Increases from zero, I.e. from the 
two-dimensional limit. Figure 6.2 shows that decay rate passes through 
a minimum, and therefore that the slowest decaying (fastest growing) 
disturbance 1s at a finite wavenumber. Furthermore, In the range of 
wavenumbers where the decay rate 1s negative, as denoted. by the shaded 
region of Figure 6.2, there are disturbances to which the flow Is 
unstable. I.e. there are disturbances which grow with time. Therefore, 
this How Is unstable. The wavenumber of the fastest growing distur- 
bance, In this case about 0.54, approximates that of a finite amplitude 
disturbance whenever wavenumber 1s a weak function of disturbance ampli- 
tude, an assumption for which there Is no a priori justification but 
which 1s conmonly made nonetheless. Tr° dominant decay rates in Figure 
6.2, and also In the figures to follow, are real-valued, and thereby 
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signify exponential decay or growth of the associated disturbance. 

The influence of the flow parameters on stability is readily exam- 
ined. Figure 6.3 shows the effect of changing the capillary number. 
When the capillary number is decreased from 0.50 to 0.45, a relative 
change of 1W, the flow becomes stable to disturbances of every wave- 
number, as demonstrated by the upward translation {accompanied by some 
deformation) of the decay rate curve. Increasing capillary number by 
10% has the opposite effect; the base flow fcacomes stable over the whole 
range of wavenumbers shown In. Figure 6.3. This clearly demonstrates that 
capillary forces arc stabilizing and viscous forces are destabilizing. 
The free surface profiles of the base flow at the same three values of 
capillary number, i.e. 0.45, 0.50, and 0.55. are also shown in Figure 
6.3. The contact line and consequently the meniscus too move upstream 
as capillary nusrber increases. However, the meniscus position is much 
less sensitive to capillary number than is the decay rate. 

It Is clear from Figure 6.3 that at some capillary number between 
0.45 and 0.5 a decay rate curve would just touch the zero decay rate 
axis. That value of capillary number is £he maximum at which a stable 
two-dimensional flow could exist, all other parameters being held 
constant. 

Flow rate also affects strongly the stability of slot and knife 
coating, «i shown 1n Figure 6.4. Reducing the flow rate for 0.5 by 10% 
leads to a dramatic destabiUzation of the base flow. Increasing flow 
rate has the opposite effect. Comparison of Figures 6.3 and 6.4 shows 
that changes in flow rate affect the meniscus position and curvature 
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more than does capillary number, and this may explain the Influence of 
flow rate on stability. One result of the dependence of meniscus posi- 
tion on flow rate Is that no stagnation line appears on the meniscus due 
to a recirculation zone near the upper solid die before the contact line 
reaches the parallel section of the channel. For this reason the Influ- 
ence of a stagnation line on stability behavior is not considered. 

Figure G.5 shows that Inertia plays a stabilizing role in slot and 
knife coating. Decreasing Reynolds number from 0.2 to 0.02. an order of 
magnitude, does little to affect either the meniscus profile or the decay 
rate curve, which shows that Inertia is inconsequential in this range. 
However. Increasing Reynolds number by an order of magnitude stabilizes 
the now significantly. This demonstrates that ribbing is indeed a vis- 
cocapillary Instability. Inertia plays a role only when It is at least 
comparable with viscous forces, and then It Is a stabilizing force. 

The effects of channel divergence angle and contact angle on menis- 
cus profile and stability are demonstrated respectively in Flguros 6.6 
and 6.7. The meniscus profile changes little, aside from being trans- 
lated upstream or downstream, as channel divergence angle is varied, 
while varying contact angle markedly affects the meniscus profile, tow- 
ever, both parameters influence the How stability; increasing channel 
divergence angle and decreasing contact angle have stabilizing effects. 
6.4. Mechanism of Instability 

A simple and somewhat heuristic rationale of ribbing Instability 
«as proposed by Pitts and Greiller (1961) and reiterated by Savage 
(1976). They derived a criterion which is based on a force balance at 
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the meniscus and on the ideas of lubrication theory. However, they 
overlooked tm important considerations: the effects of nor«l viscous 
stress at the discus and the three-dimensional nature of a rlbb'ng 
disturbance. To understand better the importance of these Issues and 
ho. they relate to the mechanism of the ribbing Instability, it „ use- 
ful to consider the complete force balance at the meniscus, which I- - 
normal stress balance as recorded in Ration (5.2.34). The terns that 

enter that equation can be clas«<n«H » k 

oe classified into three categories: (1) those 

which reoresent the chanoe in th- -« 

wange m the action of pressure on the menisci's due 

to a disturbance; (t, those *,* present the change in the action of 
normal viscous stress on the meniscus due to a disturbance; and (3, those 
which .present the change in the action of capillary force on the dis- 
cus due to a disturbance. In lubrication theory the first and third of 
these contributions can be only crudely approximated, and the second Has 
t= be omitted entirely. Here, each of the contributions » the normal 
stress balance can be evaluated from the finite element basis function 

expansions of the base flow and the eigenvector associated with the 

dominant mode of instability. 

The changes In the contributions to the no™, stress balance due 
to the onset of a disturbance are computed he* by integrating the 
weighted normal stress equation, I.e. 

3 "n - 

over a portion of the meniscus j*t downst™ of the contact line, the 
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portion where most of the action occurs (see Figure 6.1). Here u rep- 

n 

resents the perturbed velocity component normal to the meniscus, n rep- 
resents distance normal to the meniscus, "p 1s the perturbed pressure, 
Ca z uU/a 1s capillary number, u being dynamic viscosity, U being web 
speed, and <j being surface tension, h" Is the perturbed mean curvature, 
N is wavenumber and C Is an Isoparametric coordinate (see Figure 5.3). 
The normal tractions are evaluated by Integrating each of the weighted 
normal stress contributions from the contact line to a point lying 0.2 
gap thicknesses downstream, which Is the distance to the node that di- 
vides the second and third free surface elements. The weighting function 
4(5) in Equation (6.3.1) Is chosen to be the sum of two Hermlte cubic 
basis functions, two of the same basis functions which Served as weight- 
ing functions In the nonaal stress equation (5.2.34), so that « Is unity 
over the first element downstream of the contact Tine and goes to zero 
at the downstream edge of the second «ree surface clement. The magni- 
tudes of the normal stress terms are set by choosing the maximum film- 
wise displacement of the contact line due to the disturbance to be unity. 
The sign preceding the normal viscous stress term 1n (6.3.1) Is positive 
sc that when 0 n is positive, stress Is In the direction of the Inward 
pointing normal to the meniscus; on the other hand, the sign preceding 
the pressure and capillary terms in (6.3.1) are negative so that when p 
or Ca'^H are positive, stress Is In the direction of the outward point- 
ing normal to the meniscus. 

The normal stress components of the baseflow Integrated over the 
first two free surface elements are shown In Table 6.1. The curvature 
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Table 6.1. Normal tractions on the meniscus r*ar th* ******* n 

the slot or knife coat.ng nS*™ fnXTJ.u!" 6 ° f 



Traction 


Contribution 


Capillary (outward) 


0.628 


Pressure (outward) 


- 0.383 


Normal Viscous (inward) 


0.245 



Traction 


Contribution 


Capillary (outward) 


- 0.779 


Pressure (outward) 


- 0.597 



Normal Yiscous (inward) | 



1.477 
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is positive so that the net capillary force Is acting 1n the direction 
of the outward pointing normal to the free surface. The net forces due 
to pressure and viscous stress acts in the opposite direction and coun- 
terbalance the capillary force. It 1s clear from Table 6.1 that normal 
viscous stress near the contact line In this flow Is of the same order 
of magnitude as pressure, a feature that Is missing in lubrication 
theory. 

Table 6.2 shows how these meniscus tractions are changed by the 
dominant two-dimensional disturbance mode when the meniscus moves down- 
stream. The changes are per unit displacement of the contact line in 
the f1lmw1se direction. Each of the normal tractions is decreased by 
the disturbance and this leads to the following consequences: (1) the 
capillary traction on the i^niscus In the direction of the outward 
pointing normal to the free surface is decreased; (2) the traction 
exerted by pressure in the direction of the inward pointing normal to 
the free surface is Increased; (3) the normal viscous traction, acting 
in the same direction as pressure, Is decreased. These three changes 
occur In concert so that the overall balance of forces on the meniscus 
Is maintained, as It must be. Table 6.2 shows that the change in normal 
viscous stress along the meniscus due to the dominant two-dimensional 
disturbance mode is of the same order of magnitude as the change in pres- 
sure, a feature that is not included In lubrication analysis of stabil- 
ity. The realtlonshlp between the changes in the normal stress compo- 
nents and the stability of the flow is discussed further 1n the following 
paragraphs. 
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Figures 6.8 and 6.9 show how the changes In the components of normal 
traction due to the dominant two-dimensional disturbance mode depends on 
the flow parameters (excluding Reynolds number). The ordinate values are 
differences from the values given In Table 6.2. Corparison of Figures 
6.8 and 6.9 with the growth rate curves contained In Figures 6.3-6.7 
shows that the changes in the three normal tractions always turn «on> 
negative as the flow becon.es less stable. (The trend as Reynolds number 
changes agrees with this observation as well.) Thus, two-dimensional 
disturbances are least stable when they reduce capillary force (outward), 
Increase pressure force (inward), and decrease normal viscous force 
(Inward, i.e. normal viscous force is tensile). How these trends account 
for the los* of stability of the base flow to a two-dimensional distur- 
bance can be explained as follows. For the contact line to move down- 
stream along the solid die. liquid mst be driven from downstream of the 
contact line to Innately upstream of it. as shown In Figure 6.1. The 
".ore negative pressure becomes near the contact line as It moves down- 
stream, the greater is the driving force for liquid to be deliver** 
there. Normal viscous force plays a role as well. The more negative 
the normal viscous force (i.e. the less tensile) on the meniscus, the 
greater the driving force for liquid to be delivered there. This last 
point can be understood by considering an extenslonal now such as an 
accelerating liquid curtain. At the mldplane of the curtain the normal 
viscous stress along the mldplane Is tensile, consequently the normal 
viscous stress perpendicular to the raidplane Is compressive. The corn- 
Passive (negative) viscous stress normal to the mldplane corresponds to 
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liquid flowing toward It. The same 1s true at the meniscus In slot and 
kn^fe coating. Thus, two-dimensional disturbances are least stable when 
they decrease pressure and normal viscous stress at the meniscus as the 
contact line moves downstream. 

An Identification of the mechanism of instability 1n slot and knife 
coating must also account for the fact that the fastest growing distur- 
bance mode is always three-dimensional, as shown In Figures 6.3-6.7, 
Three-dimensional disturbances *re characterized by waverumber jreatar 
than zero, and are represented by Fourier basis functions in the trans- 
verse coordinate (see Chapter 5). Therefore, three-dimensional distur- 
bances of a given wavenumber u are transversely periodic and have wave- 
length 2n/N. Whatever the changes In normal stress components due to a 
three-dimensional disturbance at a plane of specified transverse coordi- 
nate, the changes are of opposite sign at a parallel plane separated by 
a distance of a half-wavelength. Likewise, whatever the change in menis- 
cus profile at a plane of specified transverse coordinate, the change Is 
of opposite sign at a parallel plane separated by a half-wavelength. It 
is enough therefore to consider how a three-dimensional disturbance 
changes the base flow 1n one plane perpendicular to the transverse coor- 
dinate axis, say at a plane that intersects the contact line where it 
protrudes downstream the furthest, for the changes at every other plane 
parallel to the first could be construct form the Fourier representa- 
tion. 

Figure 6.10 shows how changes 1n normal stress forces on the menis- 
cus due to the dominant disturbance mode depend on wavenumber for the 
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flow illustrated in Figure The ordinate values for each of the 

tractions are again differences from the corresponding values In Table 
6.2, which are for the two-dimensional limit. The following trends can 
be observed as wavenumber Increases: (1) the change 1n pressure force 
on the meniscus due to the dominant disturbance mode first decreases 
slowly, then Increases rapidly; (2) the change In normal viscous force 
increases monctonlcally; and (3) the change In capillary force - which 
equals the change In total force on the meniscus. I.e. viscous plus 
pressure force - reaches a saximum then decreases *5a1n. The maximum is 
of tensile force on the meniscus and ccrrc-porvls approximately to the 
same value of wavenumber as does the fastest-growing disturbance mode. 
The a&ove trends hold at transverse coordinate planes where the contact 
line protrudes furthest downstream; the opposite trends hold at trans- 
verse coordinate planes where the contact line Invades furthest upstream. 
Thus, the fastest growing ulsturbance occurs at the wavenumber which 
maximizes the difference in tensile force on the meniscus between planes 
where it protrudes furthest and planes where 1t invades furthest. This 
difference in stress fuels the growth of the three-dimensional distur- 
bance because It acts as the driving force for a transverse flow from 
the planes where the meniscus Invades to the planes where it protrudes, 
and thereby causes the disturbance to grow. 

6.5. Oiscussion 

The results of Section 6.3 show that the method of employing 
GaleHcin/ finite element analysis coupled slth Fourier normal mode analy- 
sis described In Chapter 5 1s a practicable means to determine the stabil- 
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tty of a steady, two-dimensional, free surface flow to three-dimensional 
disturbances. Furthermore, the eigenvectors w*1ch result from the 
Inverse subspace Iteration strategy described In Chapter 5 and Appendix 
A contain valuable Infection about the d^nant disturbance aodes. as 
shown in section 6.4. The resets of Sections 6.3 and 6.4 ^ „„ « 
*hy no. parameters Influence stability, and a,so why the fastest growing 
disturbance Is three-dimensional. 

The decay rates in Figures 6.3-6.7 show that the ribbing Instability 
<n slot and knife coating „ , nduced by ^ (l) ^ 

capillary number; (,, deceasing How rate; (3, increasing contact angle, 

-*W the solid die less lettable to the ll q uid; W decreasfng 
channel divergence ang,e; and (5) decreasing Reynolds number. Figures 
6.8 and 6.9 show that when parameters are changed In the above directions 
the changes in no™, tract1ons Kar ^ ^ ^ ^ ^ ^ ^ ' 

tnant two-dimensional disturbance .ode all become more negative - , e 
the capillary force ,„ ^ ^ . $ ^ ^ 

ecus force in the inward direction Is decreased, and pressure force in 
the inward direction is increased. The decreases 1„ no^T tract1ons 
«.ake a recirculate* How pattern such as the one shown in Figure 6.1 
"ore likely to grow, as explained in Section 6.4. and thus lead to 
Instability. 

The results ,n Figure 6.10 explain why the fastest growing (slowest 
decaying, disturbance is at a finite wavenumber. Hnen wavenumber has 
approximately the same value associated with the fastest-growing distur- 
bance mode, the net tensile force normal to the meniscus, i.e. 28u /a„- p 
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Is a maximum at transverse coordinate planes where the contact line pro- 
trudes furthest downstream and a minimum at transverse coordinate planes 
where the contact Hne Invades furthest upstream. To understand the 
connection oetween the ten^ie fcrce at the meniscus and the growth rate 
of a disturbance, 1t Is useful to examine the component of the linearized 
momentum equation (5,2.4) that Is normal to the meniscus. 

, 3u n 3u « a aui , au 1 au 1 

L 3t * n J an 1 an p J a* l a* 7T J 

Here n and i refer respectively to distance normal and tangential to the 
undisturbed meniscus, and primes denote a disturbance to the baseflow. 
Only one convectlve term appears In (6.5.1) because the baseflow has no 
velocity component normal to the meniscus. The last two terms 1n (6.5.1) 
disappear on the meniscus because of the shear stress boundary condition 
(5.2.22). Furthermore, the second term on the right side must be small 
near rhe contact Hne, for othervrfse the contact angle would be per- 
turbed by the disturbance. Therefore, on the meniscus near the contact 
line (6.5.1) reduces approximately to: 

The right hand term Is due to the exponential time dependence of the 
disturbance, i.e. e* Xt . It follows that when the tensile stress on the 
meniscus, I.e. 2 3u'/3n - p' there, 1s more positive at locations where 
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the decay rate curves fn Figures 6.3-6.7. and meniscus curvature at the 
contact line of the base Hew change as a parameter increases. Their 
hypothesis is in agreement with the trends as capillary number, How 
rate, and Reynolds number are varied. However, their hypothesis is 
contradicted by the trends as channel divergence angle decreases (B 
increases) and as contact angle increases. In the case of channel diver, 
oence angle the flow becomes less stable while the neniscus curvature 
remains almost constant; fn the case of contact angle the flow becomes 
less stable as the meniscus curvature decreases. 

!t is notable that the capillary number at onset of instability as 
predicted by finite element stability analysis (see Figure 6.3). which 
is approximately 0.5. is close to the value of 0.45 where ribbing was 
first observed in the visualization experiment with the .cylindrical 
spreader pictured in Figure 1.7. Furthermore, the wavenumber «asured 
from Figure 1.7c is about 0.53. which is only » less than that pre- 
dicted in Figure 6.2. it should be noted, however, that the configura- 
tion of the gap in the experiment was similar but not identical to that 
Illustrated in Figure 6.1. and also that How rate and contact angle 
-ere not measured in the experiment. Kevertheless, the close agreement 
between the results of the stability analysis and the experimental obser- 
vations appears to confirm the accuracy of the finite element/Fourier 
approximation. 

Figure 6.12 compares the decay rate results predicted by lubrica- 
tion theory as described In Chapter 3 with those of finite element sta- 
bility analysis. The lubrication analysis Is for flow rate specified at 
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the Inflow bounds, which is taken to be asymptotically far upstream, 
and for Coyne-Elrod type conditions at the separating meniscus. The ' 
meniscus conditions had to be computed from the finite dement base Ho, 
solutions because the contact angle, which Is 125.3', is nuch larger than 
those tabulated by Coyne and Elrod. The meniscus conditions that were 
used are given in Table 6.3. Lubrication theory does not predict the 
onset of instability at a capillary number close to 0.5. but Instead 
predicts that the How Is stable. It would be interesting to see at what 
capillary number lubrication stability analysis does predict onset of 
instability. However, 1t was not possible to find out because the con- 
tact line In the base How generated by Galericin/flnlte elements reaches 
the parallel portion of the channel at capillary numbers only slightly 
higher than 0.55 and therefore it was impossible to generate the appro- 
priate nsnlscus conditions to use in the lubrication analysis. The 
reason for the failure of lubrication theory to correctly predict the 
onset of instability i s „ wst certa1nly ttat ft ^ ^ ^ 

normal viscous stresses or, the discus, which are shown to be important 

In Tables 6.1 and 6.2 and also in Figures 6.8-6.10. 

Finite element stability calculations were also made for the expo- 
nential slot coater for the case where the contact line Is pinned and 
for capillary number ranging from 0.01 to 1. Reynolds number ranging 
from 2 to 20, and flow rate ranging from 0.35 to 0.5, but no loss of 
stability was uncovered. It is notable that at capillary numbers near 
unity and at low Reynolds number the dominant disturbance mode is pre- 
dicted to be time oscillatory when the contact line is pinned, in con- 
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trast to the situation when the contact line Is free to move along the 
die. If this mode were excited, the result would be cross-web barring. 
Hlgglns (1980) found by lubrication analysis that the extrusion coatcr, 
of which the slot coater 1s the downstream portion as shown 1n Figure 
1.1, can suffer temporal oscillations which lead to cross-web bars. 
Thus the stability of slot and knife coating appears to be sensitive tc 
the physics that governs the contact line. The best avenue for verify- 
ing hypotheses of contact line physics might be to compare experimental 
data with predicted stability behavior rather than with predicted base 
flow behavior. 

Table 6.3. Coyne-El rod type boundary conditions used to generate the 



growth rate curves shown 1n Figure 0. 12. Boundary conditions 
were computed from finite element calculations for a Reynolds 
number of 0.2, flow rate of 0.5, contact angle of 125.3° , 
and capillary numbers of 0.45, 0.5, and 0.55. 



H 



ca 



d/h. 



R/d 



0.45 



2.69 



0.821 



0.5 



2.39 



0.791 



0.55 



2.15 



0.763 
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POSTLUDE 

Mow for the first H«e, questions a5oUt ^ stsMV ^ cf a ^ 
tU^y compacted coating flow -an exanple 0 f a viscous free surface 
flc-can be answered with confidence. Ho longer is It necessary to 
Ht trust in the lubrication approximation for lack of a reliable 
approach. The key „ a f1nlte e1ement/Fo(Jrler of posslbie 

aodes of instability. Finite ele ren t analysis opens the way not only to 
determining the stability of a steady, two-dlmenslcnal base flow, but 
also to understanding the aechanls, of Instability. Such understanding 
-111 -to ft possible to avoid, or at least control , perfections and 
^iplent fll- breakdowns whenever they a* not absolutely Inherent In a 
coating system. 

The methods of finite element stability analysis developed here are 
novel. They are applicable to the entire class of steady two-dimensional 
"ows with or without free surfaces. While the development here 1s for 
Newtonian liquids with unlfonn parties, the methods can be generalized 
f ' • <• of non-Newtonlan i iquMs and * nows ^ ^ 
composition are nonuniform The feasibllty and usefulness of finite ele- 
«ent stability analysis 1s demonstrated here for s.ot and knife «.«„. 
This work is the first of Its kind, although it Is paralleled by ^ " 
ongoing research of R Uschalc {19a2) 0 „ ^ ^ ^ ^ ^ ^ 

(1982) on partially confined flows. Their analyses, however, are 1-ss 
general In applicability, as discussed below. 

Perhaps the biggest uncertainty of the demonstration of finite ele- 
-* stability analysis ,n slot and knife coating is the o,est1on of how 



216 



fine a tessellation ts needed to accurately represent the dominant dis- 
turbance modes. Finer tessellations than the ones used here were not 
possible, owing to the limited memory of the Cybers 74 and 730. The 
effects of the tessellation were explored instead by varying the size of 
the forming zone and by refining the tessellation 1n the vicinity of the 
contact line. The result: Indicated that the tessellation used here was 
fine enough to limit error In the estimates of the decay rate of the 
dominant mode to less than 5?. However, the effects of tessellation 
refinement should be tested further. This testing Is now made feasible 
at the University of Minnesota by the acquisition of a vector-processor, 
the CRAY-1. 

Finite element stability analysis can be applied to each of the 
flows illustrated In Figure 1.1. The first extension of this thesis is 
by the author's colleague Dennis Coyle, who 1s beginning to Investigate 
the stability of roll coating and has plans to investigate the stability 
of reverse roll coating as well. Coyle will be carrying out most of the 
computations o« the CRAY-1, so that he will be afforded the opportunity 
to further test the effects of tessellation refinement. Another exten- 
sion 1s being planned to slide coating, which 1s one of the most widely 
used methods In the manufacture of photographic films. 

The uses of finite element stability analysis of coating flows are 
many. It Is possible U predict just what influence changes In each of 
the operating parameters and liquid properties will have on stability. 
Complete parametric studies can be carried out to Identify stable oper- 
ating regions or "windows 0 . Furthermore, change? in die design, such as 
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the channel divergence angle in the .i.t or knife coater, can be inves- 
ted without recourse to pilot plant expectation. Indeed, new 
methods to apply coatings can be serened first by numerical analysis, 
so that less experimentation is needed and fewer dies need be machined. 

Finite element stability analysis can also be used to probe the 
mechanisms of instabilities wherever they occur. For instance, the 
importance of non*l viscous stress on the interface to the mechanism 
*f ribbing in slot and knife coating is shown here for the first time. 
Previous analyses using lubrication theory have ignored normal viscous 
stress altogether. Another important example is in slide coating. Two 
postulates have been proposed concerning the danism of ribbing in 
slide coating: one is that instabii* is viscocapillary in origin, and 
so is analogous to ribbing in slot or knife coating; the other is that 
lability is akin to the centrifugal instability in film Hows uncov- 
ened by Schweizer and Scriven (1982). Evidence was recently presented 
by Saito et al. (1982, that the mechanism is viscocapillary in nature, 
but the. issue is not completely settled. The controversy could be 
resolved once and for all by a finite element stability analysis of 
slide coating. 

Finite element analyses of steady base flows and their stability 
would not be possible without the modern high-speed digital computer. 
Without it. the sophisticated matrix calculations needed to solve effi- 
dently the large systems of algebraic equations to which base flow and 
stability analyses lead would be completely impracticable. Stability 
calculations are carried out efficiently on computers such as the 
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University of Minnesota Cyber 730 which can rapidly access at lca>l 
200,000 words of memory. For example, one stability point of the slot 
and knife coating flows reported here usually took about one minute of 
central processor time.. Virtual memory machines, such as the VAX and 
same of the newer IBM machines, could also be used, but would normally 
require longer computation times. Vector machines, such as the CRAY-1, 
are probablly ideal for stability analyses because of their large mem- 
ories, high speeds, and vector capabilities. Steve Klstler and Dennis 
Coyle are already exploiting the University of Minnesota CRAY-1 in their 
steady two-dimensional flow calculations of curtain coating and roll 
coating, respectively, and Haroon Kheshgl is doing so In unsteady two- 
dimensional film flow calculations. However, even with the CRAY-1 sta- 
bility calculations will remain relatively expensive. Therefore it Is 
imperative, especially for university research, to develop efficient 
strategies and algorithms. 

A useful strategy to treat more efficiently most flow problems that 
are defined on very large domains is to divide the entire domain up into 
those zones, if ar\y, in which the How Is nearly rectilinear, or other- 
vise close to an exact solution, and those In which It Is not. The idea 
fs to describe the nearly rectilinear Hows by asymptotic analysis and 
to describe the more complicated flows In the remaining zones by finite 
elenient analysis. The asymptotic equations are often linear, and when 
they are they can be solved once and for all. Matching conditions 
between asymptotic zones and zones In which finite element analysis is 
needed lead to vector Robin conditions which can be imposed as natural 
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boundary conditions In the finite element foliation. The result Is 
that the finite element doaain can be shrunk from what it would be. were 
".ore conventional Oirichlet or Neumann conditions used at the open-flow 
boundaries, and s, the number of algebraic equations and unknowns can be 
reduced. I„ other words, the power of the finite elenent method can be 
concentrated where it is roost needed. The final outcome is that com- 
putational costs can be cut without sacrificing accuracy. This novel 
strategy is used in Chapter 4 to represent two-dimensional base Hows 
and in Chapter 6. to represent unsteady three-dimensional disturbances. 
The same strategy should prove useful in analyzing rarv otne r free sur- 
face How problems, and in particular those In Figure 1. 1. Kistler 
(1982) has developed an analogous method for Jets and filing curtail, 
except that the asymptotic analysis there is nonlinear so that it has to 
be carried out simultaneously with the finite element analysis of the 
zone where the curtain is formed. 

Calculating stability by the approach followed here requires solving 
a large matrix eigenproblem that is asymmetric, as brought out in the' 
Prelude. To the author's knowledge, this thesis puts Stewart's 1978 
algorithm for asymmetric eigenproblems to practical use for the first 
time. Stewart's algorithm Is a relatively sophisticated method to track 
the most dangerous .disturbances in a stability analysis. Used In con- 
cert with a bil^ar mapping, a simple transformation borrowed from the 
■nethods of complex algebra. Stewart's algorithm can be used to track 
time oscillatory modes, even if they have relatively high fluency (1„ 
which case they are harder to find). Stewart's algorithm was implemented 
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here by using continuation, i.e. by specifying a previously determined 
set of dominant modes as Initial guess, to improve efficiency. The 
elgenanalysls usually required about 2/3 of the computer time needed for 
one stability calculation, or about 40 seconds o.^ \he University of 
Minnesota Cyber 730, 

Although the lubrication approximation has been used by others for 
stability analysis (see Pearson I960, Pitts and Grelller 1961, and Savage 
1977; sse also Taylor 1963) the Influence of boundary conditions on sta- 
bility predictions Is established here for the first time b> iieans of a 
side-by-side comparison. The results (see Figures 1.2 and 1.8) show that 
the Influence is strong Indeed. Furthermore, the Coyne-Elrod boundary 
condition, which appears to be the most reasonable of those compared, 
disagrees quantitatively with the finite element stability analysis (see 
Figure 1.9). Consequently, lubrication theory cannot be relied on for 
stability predictions. Moreover, lubrication theory 1s limited to flows 
1n which the domain Is nearly a parallel channel, and thus ft would be 
completely Inadequate for flows such as slide coating. 

Several possibilities exist for Improving or extending the stability 
analysis presented here, and some of them are In the offing. One Im- 
provement is to discard the independent free surface representation used 
here 1n favor of one that is consistent with the Isoparametric mapping as 
devised by Kistler (see Mstler and Scriven 19S2), This Idea is already 
being applied to finite element stability analysis by Dennis Coyle. By 
exploiting the CRAY-1, he Intends to examine more thoroughly than was 
possible in this first work the sensitivity of the dominant modes and 
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algorithm to evaluate the stability of a free surface How. 

One major aspect of a full bifurcation analysis, I.e. branching of 
families of flow states. Is missing from this work. Because the bifur- 
cating solution families are three-dimensional, a three-dimensional 
finite element algorithm would be needed to track those families. The 
present analysis is only for small -amplitude, three-dimensional distur- 
bances, which can be represented In the transverse direction by Fourier 
normal modes. Either standard three-dimensional u br1eic° elements could 
be ustd, or perhaps special ones could be constructed to represent the 
ribs with fewer degrees of freedom. In either case, the finite element 
domain would need be of only half -wavelength's width in the transverse 
direction. The possibllty of solving three-dimensional flow problems at 
reasonable cost is Just coming into view with the advent of larger and 
faster computers. Thus it may soon be possible to trace out the whole 
bifurcation diagram for coating flow problems. 

The main achievement of this thesis Is to originate and advance 
methods to evaluate the stability of coating flows, and more generally 
viscous flows of all kinds, by finite element analysis. Finite element 
stability analyses should prove to be powerful tools in designing, opti- 
mizing, and controlling coating devices, as well as understanding better 
the ways that flows lose stability. 
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APPENDIX A 

Stewart's Algorithm for th e Asymmetric Hg enjiwihi— 

In two papers Stewart (1976 and 1978) presented Ms algorithm for 
finding the m-dlmenslonal dominant, invariant subspace of a real, asym- 
metric. M x N matrix, a -or equivalent!,, finding the suhspace spanned 
by the m dominant eigenvectors. 2 { . which satisfy 

-^til iml '* - (A.l) 

where 1^1 > 1^1 > ... > U||| . In tte „ nt papep hfi ^ ^ ^ 
n.athenatica! foundations; 1n the second, the practical details. Here. 
Stewart's algorithm is described, with emphasis on the main Ideas and' 
algorithmic structure. 

Schar's theorem (see Stewart 1973 for proof), the cornerstone upon 
which Stewart's method is built, states that for every square complex 
matrix A there esists a unitary matrix. X\ the columns of which, xj, 
£ are fr*" 1 " wtors of A. i.e. vectors such that 

r-i'ts • r-iM-c" (A . 2) 

is an upper-triangular matrix. The superscript H indicates the Hermltian 
transpose (i.e.. the complex conjugate of the transpose) and C NxN signi- 
fies the space of N x N .-omplex matrices. Obviously Equation (A.2) wtt 
hold under the restriction here that A be real. The un1ta,y matrix. X'. 
may additionally be chosen so that the* diagonal elements of S\ whldT ' 
are the eigenvalues of A. appear In descending order of absoTute value 
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(I.e., absolute value of the square root of the product of an eigenvalue 
with Its complex conjugate). This ordering is assise-* in the following 
discussion. A method for finding the ordered trt* . ix S* and 

Schur matrix X' Is given below. 

The Schur vectors of a non-Hermltian matrix (1 -wtrix for 

which the Hermitlan transpose Is not equal to the matrix; enjoy many 
of the properties of the eigenvectors of a Hermitlan matrix. They form 
an orthononmal basis set and are unique provided that the absolute 
values of the eigenvalues of A are all distinct. If they are not and 
Uj^l > 1^1 - 3 ... = Ujl > U J+1 I, then the Schur vectors 

x r x j + p**-» x j are not unique, but the subspace spanned by them Is. 
Provided that |x f | > U^l, the first r Schuv vectors span the r- 
dimcndonal dominant subspace of A, I.e., the subspace spanned by the 
eigenvectors corresponding to the r eigenvalues which are largest In 
absolute value. This is easily shown by multiplying Equation (A.2) by 
X' to get 

AX' * X ' S ' , (A.3) 

for If XJ, = (x{,x£,...,xj.) and if SJ, is the r x r leading principal sub- 
matrix of the upper-triangular matrix S 1 , 1t follows that 

mm 

A^'S-S; . (A.4) 

which validates the above assertion. In the case where the absolute 
values of the eigenvalues of A are distinct, it follows from the triangu- 
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larlty of £ and the arbitrariness of r that: (1) the first Schur vec- 
tor 1s the first eigenvector of A; (2) the second Schur vector Is the 
linear combination of the first two eigenvectors of A that makes It 
orthogonal to the first Schur vector; and (3) the k-7h Schur vector Is 
likewise that linear cognation of the first k eigenvectors of A that 
makes it orthogonal to the first k-1 Schur vectors. Because of tte 
uniqueness properties of the Schur vectors, the third point holds for 
X' whenever U^l > ^ > ^ eve „ ff ^ ^ ^ ^ 

eigenvalues of A are not all distinct. 

The elgenproblero is practically solved once the solution to equa- 
tion (A.5) is found, because the r dominant eigenvalues appear along the 
diagonal of s;, and the corresponding eigenvectors can easily be com- 
puted by solving the r x r eigenproblem 

Ir^'H^k k"l,2 r (A . 5 , 

fery^ and then premultlplylng by to get 

* irA • k - 1.2 r (a.5, 

Equations (A.4) and fA.5) imply Equation (A.6) because 

2£&-V2W\VH). (A.7J 

so that {£ j£) ast be the eigenvector of A corresponding to the eigen- 
value A k> and therefore It mist be the same"as z fc . 

The original problem of finding the dominant subspace of A is now 



reduced to tao subproblems. First, Equation (A. 4) roust be solved to 
find the Schur matrix, V , and the upper triangular matrix, S 1 . Then, 
Equations and (A.6) have to solved tc find the eigenvectors, z . 
A difficulty with this strategy 1s that complex arithmetic has to be 
used extensively whenever two or more of the eigenvalues are complex 
conjugate pairs. Furthermore, solving Equation (A.S) requires using an 
algorithm to find the eigenvectors of a conplex matrix. As a result, 
the computer memory required to store X' and S' would be doubled and the 
computation time needed to find the dominant subspace of A would be 
significantly greater than in the situation where only real arithmetic 
Is required. Fortunately, complex arithmetic can be avoided entirely by 
working with what Stewart calls the "quasl-Schur fora." The triangular 
• 3^. is traded for a quasi-triangular S^ 

5r"£iir • ST 6 *™ • L- (A - 8) 

and the matrix of quasi -Schur vectors, = (x^x^,-..)^), needs only to 
be orthogonal Instead of unitary. By quasi-triangular Is meant that 
either I x I or 2 x 2 block matrices appear along the diagonal; the I x I 
and 2 x 2 blocks contain real eigenvalues and ccsiplex conjugate pairs, 
respectively, Tne existence of a solution to Equation (A.S) for arbi- 
trary real A was proved by Wilkinson (1965). Schur 1 s theorem can be 
applied tc ^ tc show that 
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and therefore that 
X* « x X* 



(A.10) 



The first equality In Ration (A.9) follows from t>: uniqueness proper- 
ties of the triangular matrix Equation (A.IO) establishes that the 
Schur matrix, i s Just ^ product of ^ orthogOMj quas1 . Schur 
aatrix, X., and an r x r unitary matrix, x\ The first equality 1n Equa- 
tion (A.9) also establishes that 



S . » X* X* H , 



(A. 11) 



so that the quasi-traingular £ can always be constructed by a unitary 
transformation of the triangular SJ,. 

In practice, S,. can be found directly, and thus need never be 
obtained. This Is what Stewart's method does, and thereby It avoids 
complex arithmetic. The eigenvectors of A can be found once ^ and X 
are known by using the quasi-Schur equivalents of Equations (I^S) anf 



lrh'\h k - 1.2 r 



h'lrh. k = 1,2,. ...r 



(A. 12) 
(A.13) 



• The next few paragraphs sketch out the three main steps fn Stewart's 
algorithm for solving Equation (A.9). a simple numerical example serves 
to Illustrate each of the main steps. Subseouent paragraphs describe how 
these three main steps and interconnecting logic for* the overall algo- 
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rlthm. The final paragraph summarizes some of the practical aspects of 
Implementing this ilgorithm. 

The power step 1s, 1n the Iterative procedure of the algorithm, the 
most frequently executed, and In a sense It 1s the most Important of the 
three main steps. It 1s solely responsible for the convergence from an 
Initial subspace represented by the matrix qJ 01 = ( 5 j 0) , £ Q) . ^° 5 ) 
to the dominant sui/Space, which 1s represented by the quasl-Schur matrix 

£r = £g i^* P° wer step Is amazlnoly simple: It Is 

little more than a matrix multiplication. In the case that only one 
elgenpalr Is sought, the v-th Iterate £ (v| defined by 



converges to z^ as v Increases provided that Is nonzero In the expan- 
sion 



(A. 14) 



(0) 



N 



1 




(A.15) 



and also that the first eigenvalue of A dominates: 



Ujl > |A 2 I > IjU > ... > U N I . 



(A.16) 



This 1s easily shown, because Equation (A.14) Is equivalent to 




N 



Vl + & 7 1 'VV *| 



(A.17) 



and U^/Xjl < 1 for 1 « 2,3 H. Equation (A.14) cannot be used In 
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practice because ^ 1s not known a prion. Furthermore, computation of 
A v requires an inordinate amount of work, about M 3 multiplications and 
additions per Iteration when A is full. These problems are overcome by 
replacing Equation (A.16) by 

* - 4 /c » (A.18) 

wler, e (w » is the Euclidean nonn ofA^', i.e., >> 1s ^ ^ 

of the sum of the squares of the absolute values of the e7** Ilts ,„ 
*-i " • In practice, the normalization need not be done at each itera- 
tion. For comparison, Equation (A.18) requires only about M 2 ♦ N multi- 
plications and additions per iteration. Unfortunately. ^ *n, not 
converge when the first two eigenvalues have equal absolute values. For 
example. If the eigenvalues are dominated by a complex pair so that tin 
first Inequality In Equation (A.16) does not hold, but 

does hold, then the iteration defined by 

Q^ +i) * a o (v, r c (v, i- 1 

* l - J ' (A.20) 

where ^» . (jM ^1, and ,W |$ d1agonal ^ ^ ^ 
to the Euclidean nonns of a£\ a „o a^>. converges to the donrtnant 
subspace. However, convergence of Q* w) to the dominant subspace cf A 
cannot be guaranteed even when |x r l > Ix^l if ttnlt e decimal artJ 
metlc is used. Consider the third-order matrix - 
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A = 



3 0 0 

4 Z 0 
0 0 1 



(A.21) 



which has eigenvalues and eigenvectors 



If tiie Initial guess is 



10 0 
4 1 0 
G G 1 



(A.22) 



"1 


1 " 






1 


0 


f 


(A.23) 


1 


-1 







equation (A.21) gives the result (without normalization): 



■~3 U 



of. 



4.(3 V -2 U ) ♦ 2 W 4.(3 V -2 V ) 
1 -I 



Taking the limit v*» and normalizing give: 



(A.24) 



2 \h-» 2 



i i 

4 4 

0 0 



/ m 



(A.25) 



Both vectors converge to the dominant eigenvector. Zj! This could have 
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been forseen Inasmuch as 1^ > ^ inC , oth ^(0) and ^ 
component of the first eigenvector. In finite decimal arithmetic the 
problem Is e^on more severe; after a finite number of steps a J v) and ^ 
become identlcsl. This shortcoming Is remedied by periodically ortho- 
nonaallzlng the columns of Q^>. This orthonomallzatlon is the second 
main step In Stewart's algorithm. 

A fir «- Schl| rfdt orthonormallzation <t» r U ^ ^comi step. In It 
the columns of gj. u) are .lormalized and made orthogonal: 

Q (v) ■ 0 (v) (R (v, r l 

* 'ir ' (A.26) 

where rJ v1 is and r x r upper triangular matrix chosen by a modified 
Gram-Schmidt algorithm to orthonomlli. (see Stewart 1973 and, for 
a program, Rutlshauser 1970). The second step need not be executed after 
eart power step, but just frequently enough to assure that the columns 
of 0/ do not become linearly dependent. While the Gram-Schmidt method 
Is less stable numerically than Householder's method (see Stewart 1973), 
it is computationally efficient and simple to program'. In practice, 
numerical instability is restrained by adopting a conservative criterion 
for the frequency to perform the orthonormallzatlon step. This criter- 
ion Is described below. 

In the above example an orthonormalizatlon of q| v) before taking 
the limit w*» produces the result: 



" 1 4 
4 -1 
0 0 



' /17 (A.27, 
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The column vectors of Q 9 do span the dominant subspace of A, and In fact 

1s identical to the quasl-Schur matrix Xj. The dominant eigenvalues 
are found by forming the unitary projection of the matrix A Into Q^, as 
defined 1n Equation (A. 8): 



(A.28) 



The diagonal elements of Sg are the dominant eigenvalues of A, which are 
recorded in Equation (A. 22). The eigenvectors of A are found by solving 
Equation (A.12) for (£ r yJ and then solving Equation (A.13) for (i^Zg). 



ri 4A'ifi 

Lo -1//T7J 



u vh ) • 



1//17 4//17 
4//T7 -1//17 
0 0 



p 1 4//17 T 
Lo -1//17 j 



1//17 C 
4//17 1 
0 0 



(A.29) 



(A.30) 



These are indeed the first and second eigenvectors of A, which are 
recorded in Equation (A. 22). 

The above solution was obtained as a limit, I.e. by letting the 
iteration count v increase without bound. A numerical algorithm on 
the other hand must terminate after a finite number of Iterations. 
Obviously, a convergence criterion has to be constructed. This leads 
to the third main step In Stewart's algorithm, which he calls the 
Schur-Raylelgh-Ritz step , or SRR step for short. The SRR step is exe- 
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cuted prior to each convergence test, and It can accelerate convergence 
as well. The acceleration of convergence I, not due to an ioip.-ove.nent 
In the current estimate of the dominant subspace, but rather to an 
"unscrambling" of the vectors which span the current estimate of the • 
dominant subspa«. This unscambHng of vectors amounts to Improving the 
accuracy of the more dominant quasi-Schur vector approximates by using 
Information contained In the less dominant approbates. For example, 
if an r-dlmensional subspace Is carried in the iteration process, but 
only an m-dlmenslonal dominant subspace Is sought, where «. < r, the SRR 
step reduces the number of Iterations needed to obtain a solution of 
given accuracy, as compared with the number needed without the SRR step. 
On the other hand, where m - r the number of iterations would not be 
reduced by the SRR step, but the computational work would be reduced 
nevertheless becausa of a deflation process that Is Included in Stewart's 
algorithm, as discussed below. The SRR step proceeds as follows. First, 
A is projected Into the subspace the basis of which is the approximate 
set of quasl-Schur vectors, Q^ v ': 

-* K Zr 1 Z2r ' (A.31) 
Then Schur's theorem is applied to 8^': 

T {v) » fY (v, ) T B (v > y< v > 

•* l ir > !r It (A.32) 
to obtain the set of quasi-Schur vectors, Y* v) , and 7< v) , the quasi- 
triangular unitary transformation of bJ v} . The method employed to 
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solve Equation (A. 32) 1s the well-known OR algorithm with double Impli- 
cit shifts (see Stewart 1973 and Wilkinson 1965 for a comprehensive 
description) modified to order blocks of eigenvalues as they emerge (see 
Smith et aU 1974 and Wilkinson and Relnsch 1971 for FORTRAN programs). 
The dominant algenvalues of A are approximated by the eigenvalues of the 
I x I and 2 x 2 diagonal blccks of t| v ', and the dominant quasi-Schur vec- 
tors are approximated by the columns of 

g(v) = l(v) (v) 

2r = }jf J* • (A. 33) 

Stewart (1976) showed that 5^' is a better approximation to the quasi- 
Schur matrix than Q^ v '. This fact can be understood 1n an intuitive way 
by noting that In 

(A.34) 

the columns of $£ v ' form a different basis set for the same subspace as 
is spanned by Q^ v ' ( a basis set so chosen that 1t quasl-trlangul arizes 
A. Although Equations (A.34) <?nd (*.8) look alike, they are not iden- 
tical because need not span the dominant subspace. However, as v 
increases, 9^ v) and T^ vl do approach ^ and S^., respectively. 

The benefit of Introducing the SRR step Is brought out clearly by 
returning xo the previous example. Setting v = 5 fn Equation (A. 24) 
provides an approximation to the dominant subspace of A: 
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243 243 1 
844 812 
1 -1 



(A.35) 



An application of the orthonormal izatlon stsp produces 



0.27668 0.93620 "' 
0.96096 -0.26987 
_ 0.00114 -0.21667 



(A.36) 



The accuracy of qJ- can be determined by fomlng Its product with the 
exact quasl-Schur matrix from Equation (A.27), 



(5) 



Sz'Sz 



L o. 



99937 
03535 



-0.03426 "7 
0.97564 _[ 



(A.37) 



The departure of this result from the 2 x 2 Identity matrix indicates the 
nagmtude of the remaining error. Equation (A.37) also conveys the fol- 
lowing information: (1) the first column vector In Q< 5 > is mch mre 
closely parallel to Its corresponding quasl-Schur vector than Is the 
second column of Q< 5) , as revealed by the closeness to unity of the 
diagonal elements; (2) the two column vectors of $<« contain about equal 
components of the opposite quasl-Schur vectors, as revealed by the off- 
diagonal elements; and (3) therefore the second vector In & S) contains 
a larger component of the omitted third quasl-Schur vector. The third 
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point results from the fact that the complete set of quasl-Schur vectors 
forms an orthonormal basis for the vector space. All three points are 
predicted by Stewards convergence rate analysis, which Is set out below. 

The estimate of the dominant quasl-Schur vectors 1n Equation (A. 36) 
can be Improved by taking an SRR step. First, equation (A. 31) Is solved 
fyr S^ 51 and then Equation (A.32) Is solved for J^ 5) and t£ 5} : 



R {5) 



T 3.14005 3.86636 1 

(A.38) 

L-0.03861 1.82052 J 

,c\ r 0.99948 0.03230 1 
Yi 5) ■ (A.39) 
L-O.C3230 0.99940 1 

The first column vector of Y^' Is the normalized dominant eigenvector 
of the ^cond i$ chosen to make the pair orthonormal. 

, 0 r 3.01509 3.90498 "T 
JV ] ■ - CA.40) 

L 0 i.94549 J 

The dominant eigenvalues are approximated by the diagonal elements of 
T^ 51 (cf. Equations (A. 22) and (A. 28)). Finally, t)^ 5 * Is computed from 
Equation (A. 33): 



jfS) 



0.24623 0.94665 
0.96133 -0.23869 
0.00814 -0.21652 



(A.41) 
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Preaulttplytng jj» by the exact quasl-Schur matr.x. given in Equa- 
tlon (A.27) then leads to a revealing comparison with equation (A.37) 



-* J L °-°°572 0.37o*28 J 



(A. 42) 



The off-diagonal element tn Equation (A.42) are smaller In magnitude 
than those In Equation (A.37) by an order cf magnitude. The unscram- 
bling process performed by the SRR step brings about an Improvement In 
the approximate quasl-Schur vectors by making each of then, more nearly 
orthogonal to all of the exact quasl-Schur vectors except the one to 
which It Is converging. The convergence criterion used In 
(A.37) and (A.42) cannot be used in practice because the exact quasl- 
Schur matrix Is not *,c-,n a priori. A suitable convergence criterion 
Is to require that the first m column vectors of the residual matrix 
defined by 

-J 1 J ir 1 (A.43) 
have small norms: 



^ < « far ! - 1,2 



(A.44) 



where R (v) = (>, (v > *M . (v). . 

_r *M • £4 . ^ ) and is the Chebyshev norm. 

I.e.. the maximum of the absolute values of the elements in a vector. 

Computing £ 5) for the above example gives 
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(A.45) 

The Chebyshev norms of the column vectors of R^' are 0.016*0 and 
0.17293. These Indicate the errors In the two approximate quasl-Schur 
vectors, ^j 5 * and respectively. The first Is seen to be more accu- 
rate than the second by an order of magnitude. 

Stewart (1976) proved t*»at this scheme converges provided U r l > 
U^l. Furthermore, he proved that If lx % \ > U z+1 li the subspace 
spanned by the first l vectors < n ^ V * converges to the £-41mens1ona1 
dominant subspace of A at the rate 0( |A r+1 A Jt l v ). Unfortunately, the 
requirement that U r l > U^l is necessary for convergence. In prac- 
tice, there 1s a danger that r could be so chosen that a complex conju- 
gate pair of eigenvalues straddles the cutoff. Stewart (1976) shows by 
a 3 x 3 example that this would prevent the algorithm frca converging. 
The only practical insurance against this hazard Is to increase or 
decrease the subspace dimension, r, by one when it appears that the 
Iterative procedure Is not converging. 

Figure A.l 1s a flowchart of Stewart's algorithm which shows the 
three main steps and interconnecting logic. This flowchart 1s described 
in the following paragraph with emphasis on the interconnecting Togic. 

The algorithm begins by initializing Q* 0 *, either by continuation 
or by a random process. If zeroth-order continuation Is used, the Ini- 
tial guess Is a set of vectors which span the dominant subspace at a 



"0.00372 0.03673 
R* 5J = 0.00908 0.01961 
-0.01640 -0.17293 



f. START ~) 

INITIALIZE 
0 IT » 0 



PERFORM; 
SRR STEP 



COMPUTE 




RESET I '• 



COMPUTE:. -DORT, 
NXTSRR", KXTORT 



I POWER STEP l^-? ' : j 
(IT » IT + 11 - I 



I 



UPDATE 
NXTORT . . : 
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ORTHONORM. ' 

. .STEP.. 




STOP ^ 



YES 



Figure A.I. Flowchart of .Stewart's.Algor.1«Bn. 
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nearby point in the parameter space; 1f first-order continuation 1s used, 
the Initial guess 1s generated by linear extrapolation from a nearby 
point in the parameter space, as Ascribed by Sillissn (1979! * The 
alternative to* continuation 1s to assign values to the elements of (J* 0 ' 
by using a random number generator. In either case, orthonormallzatlon 
is performed {unless unnecessary) to complete the initialization, after 
*Mch an SRR step Is taken and residuals computed. A check 1s made to 
see if the first previously unconverged column vector has converged 
(initially the first vector). If so, the parameter, L, which points to 
the first unconverged vector in£ (v> , 1s Incremented. The (L-l) con- 
verged column vectors 1n£* v * are skipped over in most subsequent opera- 
tions in order to save confutation time. This process Is called defla- 
tion (see Stewart 1976). When L > M, M being the dimension of the domi- 
nant subspace which is sought, the calculations are halted. Eigenvectors, 
1f desired, are computed 1n a post-processing step as described above. 
If L < M, the Iteration counter, I, Is compared with a user-specified 
maximum, MAX. The computations stop If I * MAX; otherwise, the parame- 
ters MXTSRR and NXTORT are calculated, these being the Iteration oumbers 
at which the next 3£ :r.d crthonormalization steps, respectively, should 
be performed. MXTSRR Is set to the minimum of three values: (1) the 
maximum number of Iterations, MAX; (2) the iteration at which the next 
quasi-Schur vector is expected to converge; and (3) twice the previous 
value of NXTSRR (set to S Initially). The iteration number on wnlch the 
next quasi-Schur vector 1s expected to converge Is comptued by extrapo- 
lation using the linear convergence rate predicted by Stewart. The 
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third criterion is necessary because linear convergence only occurs 
asymptotically. NXTORT Is set to the nrinlmum of SRR and I + DORT, where 
DORT is calculate-J f a conservative fonoula (see Stewart 1976) which 
insures that no more than eight decimals of accuracy can be lost during 
an orthonoraalizatto,, step. text, the power step is repea ted until I . 
«XT0RT. I 1s lncrenente<1 after Mch power ^ ^ orthononna1l2ation 
step is then taken, after which NXTORT 1s set to the minimum^ NXTSRR 
and I ♦ OORT. The power and orthononnalization steps are continued In 
this fashion until I- mm: The cycle Is then seated beginning 
with the SRR step and continuing until either the convergence criterion 
is satisfied or the maximum number of Iterations 1s^ reached'. 

A single call to a subroutine called SRRIT, which drives otter 
subroutines, 1s quired to activate 'Stewart's algorithm. The user « t 
specify parameters which set the number of column vectors In the matrix. 
£ . the dimension of the dominant subspace to be found, the maximum 

° f 'Z Ste " 1ter8ttons * * «... and the method of initiali- 
sation for Q ' -either by continuation br by a random assignment pro- 
cess. The writer found that the number of Iterations needed to attain" 
convergence (with W in tquatjon ( A .«„ was ^picaJly „duced by 

an order of magnitude when the initial suess MS generated by zeroth- 
order continuation rather than by random guess. A subroutine^ which 
Stewart called ATQ. must be supplied by the user to premultiply by 
A. The smallest eigenvalues (in the sense of absolute value) of" Eq uat1on 
(A.1, can be found by Stewart's algorithm by rewriting that equation as ' 
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L~ l Li " n 1 -1 ■ (A.46) 

where = Inasmuch as the largest correspond to the smallest 

A r The smallest eigenvalues of the generalized eigenprobTero fc 

ill 3 X i l±i ' U.47) 

can also be found by transforming the problem to 

I' 1 Uli 3 n f £i • (A. 48) 

Only the subroutine ATQ need be changed to solve Equation (A.48). Given 
the iterate, Q (v, t the next iterate is calculated by first premulti plying 
1t by M and then solving the linear system 

i^'if (A.49) 

for£* v+i '. Because this linear system has to be solved nary times 
during the iteration process, much computational efficiency Is gained 
by computinc the triangular decomposition of J (see Wilkinson 1965} 
before the first <terat1on. When J 1s the Jacobian matrix for a set of 
equations generated by the finite element method, it Is convenient to 
use Hood's frontal method (see Hood 1976) to trlangularlze J. Solving 
Equation (A.49) then requires only a call to Hood's subroutine RESOL, 
which needs to be modified tc accomodate multiple right-hand sides. 
The eigenvalues of Equation (A.47) found by the strategy outlined 
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above are those of smallest magnitude; the eigenvalues sought in coating 
flow stability analysis are those of anal lest real part. These two sets 
of eigenvalues - the ones found by applying Stewart's algorithm to 
Equation (A.48) and those of greatest importance to stability assessment 
-will be the same provided the eigenvalues nearest the imaginary axis 
do not have' Imaginary parts that are large compared with the magnitude 
of ether eigenvalues near the Imaginary axis (see also Chapter S). a 
bilinear mapping, described in Chapter 5, can be used to transform 
eigenvalues of smallest real part Into ones of smallest magnitude so 
that they can be Identified easily. 
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